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Summary 


Our  research  activity  has  been  set  in  the  context  of  a  nonlinear  finite  deformation  prob¬ 
lem  describing  unsteady  membrane  inflation.  We  have  developed  theory  and  prototype 
software  that  is  capable  of  estimating  the  errors  in  several  physically  relevant  engineering 
Qol's.  In  the  first  part  of  the  project  these  estimates  were  confined  to  only  the  finite  ele¬ 
ment  discretization  error.  In  the  last  stages  of  the  project  we  have  been  able  to  extend  this 
error  estimation  to  encompass  the  modelling  error  that  results  from  neglecting  the  inertia 
term  in  Newton's  second  law  of  motion.  Full  details  are  given  below;  it  is  expected  that 
these  findings  will  appear  in  one  or  more  journal  publications. 

In  this  project  we  have  applied  state-of-the-art  mathematical  theory  to  a  practical  en¬ 
gineering  problem  in  membrane  inflation.  We  have  applied  modern  methods  of  a  posteriori 
error  estimation  in  order  to  address  the  error  in  engineering  quantities  of  interest  that  arise 
from  modelling  simplifications  and  also  from  finite  element  approximation.  Our  results 
show  that  this  technique  works  in  that  it  can  be  used  to  adaptively  optimise  the  finite  ele¬ 
ment  mesh,  as  well  as  give  an  indication  of  the  physical  error  contained  in  the  modelling 
assumptions. 

The  eventual  benefit  of  research  of  this  nature  and  the  downstream  uptake  of  the  tech¬ 
niques  and  methodology  into  commercial  engineering  software  cannot  be  overstated.  By 
isolating  and  estimating  the  two  principal  sources  of  error  in  physical  simulations  (mod¬ 
elling  and  discretization)  the  software  will  allow  the  practitioner  to  assign  a  degree  of 
confidence  in  the  results  that  has  never  before  been  possible. 
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as  defined  by  the  payment  schedule  of  the  contract. 


1  Introduction 

This  project  has  been  concerned  with  the  reliable  computation  of  engineering  Quantities  of 
Interest  (Qol),  and  the  application  of  these  mathematical  methods  to  problems  of  interest 
to  the  US  Army  (informed  by  our  long-standing  contact  with  Dr  AR  Johnson,  ARL,  VTD, 
LaRC). 

Briefly,  an  engineer  is  often  confronted  with  a  difficult  problem,  the  'fine'  problem, 
which  in  practical  terms  is  an  intractable  computational  task.  Usually  it  is  not  every  as¬ 
pect  of  the  solution  that  is  of  interest  but  only  a  certain  quantity,  the  Qol.  For  example, 
this  could  be  a  localised  stress,  a  mode  shape  or  natural  frequency,  or  even  a  measure  of 
dissipated  energy  (damping). 

To  overcome  the  intractability  it  is  routine  to  replace  this  fine  problem  with  a  'coarse' 
problem  which,  it  is  judged,  contains  enough  of  the  physics  to  reliably,  but  approximately, 
capture  the  Qol.  This  entails  a  modelling  error.  Furthermore,  when  this  coarse  problem  is 
loaded  in  to  finite  element  software  the  computed  Qol  then  involves  also  a  discretization 
error. 

For  some  problems  experienced  engineering  judgement  alone  may  be  enough  to  de¬ 
termine  whether  or  not  the  computed  Qol  is  accurate  enough  for  the  purpose  of  the  ap¬ 
plication.  However,  sometimes  the  engineer's  experience  will  suggest  that  the  Qol  may 
be  computed  inaccurately,  or  subtle  effects  in  the  problem  may  give  rise  to  significant  but 
unanticipated  (and,  possibly,  undetected)  errors  in  the  Qol  due  to  either  or  both  of  the 
modelling  approximation  or  the  finite  element  approximation. 

Our  research  has  focussed  on  the  development  of  modern  techniques  in  computa¬ 
tional  mathematics  which  afford  a  computational  procedure  for  the  'black-box'  estima¬ 
tion  of  both  the  modelling  error  and  the  discretization  error  in  the  computed  Qol.  This 
computed  error  estimate  is  then  extra  information  that  the  engineer  can  make  use  of  in 
determining  whether  or  not  the  computer  model  is  adequate  at  the  level  of  need. 

The  basic  idea  behind  this  mathematical  technique  is  to  compute  the  solution  to  the 
coarse  problem  and  place  it  into  the  fine  problem's  equations.  The  extent  to  which  this 
solution  does  not  satisfy  the  fine  equations  (the  residual)  gives  a  rough  measure  of  how 
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adequate  (or  otherwise)  the  coarse  solution  is.  However,  this  is  a  non-specific  measure  of 
'adequacy'.  To  connect  this  residual  to  the  engineer's  specific  Qol  requires  the  additional 
computational  step  of  determining  a  dual  solution.  This  dual  solution  then  multiplies  the 
residual  at  every  point  in  the  finite  element  mesh  and,  by  a  prescribed  procedure,  these 
weighted  residuals  are  summed  up  and  the  result  is  an  estimate  of  the  errors  in  the  Qol. 

This  dual  solution  is  not  unlike  the  influence  (or  Green's)  function  from  applied  math¬ 
ematics.  Its  value  lies  in  the  fact  that  a  large  residual,  at  a  patch  in  the  finite  element  mesh, 
may  end  up  being  weighted  by  a  very  small  dual  solution  (over  that  patch).  The  net  result 
is  that,  although  the  solution  quality  is  low  in  that  patch,  the  contribution  to  the  error  in  the 
Qol  is  small,  or  even  negligible.  This  information  can  then  be  used  to  provide  a  guide  to 
local  mesh  and  model  refinement.  And  this,  in  turn,  allows  the  engineer  to  tune  the  model 
and  mesh  to  determine  the  Qol  accurately  with  the  minimum  amount  of  computational 
effort  on  the  coarse  problem. 

Of  course,  this  extra  functionality  comes  at  a  price  -  we  can  think  of  it  as  a  tax  on 
the  cost  of  the  coarse  problem.  The  tax  is  levied  on  the  extra  information  on  the  error 
and  is  manifest  as  the  cost  in  computing  the  dual  solution.  We  recognise  that,  in  some 
applications,  this  tax  is  too  high  to  bear  but,  in  other  applications,  it  can  be  seen  as  small 
relative  to  the  overall  cost  of  the  engineering  problem. 


2  Dissemination  activity  at  ARL,  VTD,  Langley 

This  is  a  summary  of  the  trip  to  ARL,  VTD,  LaRC  made  by  Whiteman,  Shaw  and  Warby 
to  visit  with  Dr  AR  Johnson  during  9-15  October  2007. 

The  primary  purpose  of  the  trip  was  to  disseminate  findings  from  this  ARO  research 
project,  and  also  to  effect  a  technology  transfer  though  discussion  of  our  advanced  applied 
engineering  mathematics  and  by  a  transfer  of  our  working  computer  software. 

The  secondary  purpose  of  the  trip  was  for  us  to  understand  current  needs  within 
ARL's  broader  mission,  and  to  comment  on  how  our  research  expertise  might  contribute 
toward  that  mission. 

The  next  part  of  the  report  summarises  our  activities  and  conversations  during  our 
trip  to  ARL,  and  (where  appropriate)  attempts  to  sensibly  comment  on  whether  or  not  our 
technology  may  be  able  to  make  a  practical  contribution. 

•  Seminar  presentation:  Warby  gave  a  seminar  presentation  of  our  research  on  Octo¬ 
ber  10,  2007.  The  audience  included  Dr.  AR  Johnson  as  well  as  Drs.  Bey  and  Tessler 
(NASA). 

•  Code  transfer:  Our  working  computer  software  was  transferred  to  Dr.  AR  Johnson 
on  October  13,  2007. 

•  Theory  transfer:  A  detailed  run  through  of  the  mathematics  underlying  the  Qol 
error  estimation  technique  was  made  with  Dr.  AR  Johnson  on  October  13,  2007.  It 
was  noted  that  a  comprehensive  guide  to  this  theory  will  be  contained  in  this  close 
out  report  to  ARO. 
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•  Discussions  with  Dr  RP  Thornburgh:  We  met  with  Dr  Robert  P  Thornburgh  on 
October  11,  2007.  Dr  Thornburgh  explained  the  difficulties  he  has  encountered  in 
the  reliable  computation  of  the  natural  frequencies  and  mode  shapes  (the  Qol's)  of 
helicopters  suffering  inflight  bending  loads. 

1.  He  explained  that  such  a  small  effect  as  modelling  clamped  boundary  condi¬ 
tions  as  a  continuous  line,  as  opposed  to  a  sequence  of  bolts,  can  throw  off  the 
computation  of  the  Qol  by  a  significant  amount.  Indeed,  even  variations  in 
the  diameter  over  which  the  bolts  were  clamped  introduced  a  significant  mod¬ 
elling  error.  This  had  non-trivial  implications  as  to  when  and  how  apparently 
secondary  structures  such  as  cockpit  glass  and  door  panels  should  be  included 
in  the  coarse  model. 

2.  Another  significant  difficulty  arose  in  the  modelling  of  the  tail  panels.  Mod¬ 
elling  the  panel  in  each  bay  by  a  single  finite  element  resulted  in  a  computa¬ 
tionally  feasible  problem  but  did  not  allow  the  (coarse)  finite  element  model  to 
predict  panel  buckling  (through,  for  example,  the  introduction  of  surface  im¬ 
perfections).  Adding  enough  elements  to  capture  this  effect,  however,  pushed 
the  model  into  the  computationally  unfeasible  domain. 

3.  Dr.  Thornburgh  also  introduced  us  to  the  current  interest  ARL  has  in  uncer¬ 
tainty  analysis  and  the  work  he  had  carried  out  with  LTC  Robert  B  Floersheim. 
We  will  comment  further  on  this  below  and  simply  remark  here  that  we  hope, 
in  due  course,  to  obtain  a  copy  of  LTC  Floersheim's  thesis. 

We  noted  that  while  our  technology  could  make  a  contribution  in  the  first  two  areas, 
it  is  possible  that  the  tax  would  be  too  high.  Dr.  Thornburgh's  main  requirement  is 
that  of  obtaining  an  accurate  Qol  within  reasonable  time.  The  computation  of  the 
dual  solution  would,  although  giving  estimates  of  the  error  in  the  Qol,  introduce  a 
significant  time  delay  into  the  entire  computational  procedure. 

On  the  other  hand,  it  is  conceivable  that  our  methodology  could  make  a  significant 
contribution  in  terms  of  uncertainty  analysis  (see  below). 

•  Discussions  with  Dr  A  Tessler  (NASA):  On  October  12,  2007,  we  discussed  with  Dr. 
Alex  Tessler  (NASA)  his  major  contribution  to  zig-zag  theory  in  the  stress  analysis  of 
beam  and  plate  structures.  It  is  clear  that  such  models  can  form  a  fine  problem  and 
that  our  Qol  error  estimation  is  relevant  to  this  area. 

•  Further  discussions  with  Dr  AR  Johnson:  On  October  12,  2007,  we  discussed  with 
Dr.  AR  Johnson  the  Natick  Soldier  RDEC's  ongoing  work  on  air-supported  struc¬ 
tures.  These  structures  are  susceptible  to  significant  ageing  effects,  geometric  and 
material  imperfections,  large  design  tolerances,  large  deformations,  and  can  be  mod¬ 
elled  by  various  mathematical  techniques.  In  the  context  of  uncertainty  analysis  it 
seems  that  our  technology  may  be  able  to  make  a  significant  impact  here.  For  exam¬ 
ple,  assume  that  in-service  conditions  vary  randomly  about  some  mean  (expected) 
values.  The  dual  solution  could  be  obtained  just  once  on  this  state  and  then  used  to 
determine  an  optimal  mesh  as  well  as  judge  whether  the  modelling  error  in  the  Qol 
was  at  an  acceptable  level. 

Given  this  optimal  mesh  the  air-supported  structure  could  be  simulated  efficiently 
many  thousands  of  times  with  randomly  perturbed  inputs.  Each  time  a  Qol  would 


-5- 


be  produced  and  from  these  outputs  we  could  build  its  probability  density  function. 
This  would  be  a  valuable  design  tool  in  terms  of  when  and  how  to  select  design 
tolerances  and  would  inform  quality  control  procedures  in  terms  of  variability  of 
construction.  Note  that  now  the  'tax'  is  at  the  order  of  1/1000. 

Although  this  project  is  now  at  a  close  we  are,  in  terms  of  the  above  comments  and 
observations,  keeping  in  touch  with  Dr  Johnson  to  further  discuss  the  feasibility  of  seeking 
further  ARO  funding  in  order  to  contribute  to  the  Army's  needs. 


3  The  theoretical  framework 

In  this  section  we  outline  the  basic  mathematical  framework  behind  this  error  estimation 
technique.  The  presentation  is  deliberately  'abstract'  so  as  to  avoid  giving  specific  details 
on  the  physical  properties  of  the  underlying  problem.  These  details  are  given  later. 

The  main  idea  behind  modelling-error  estimation  is  to  consider  a  'fine',  or  high  fi¬ 
delity,  model  which  captures  all  the  physics  relevant  to  the  problem,  and  a  'coarse'  ap¬ 
proximation  which,  we  suppose,  includes  only  those  physical  effects  that  are  needed  for 
the  computer  simulation.  The  issue  then  is  to  quantify  the  resulting  'modelling  error'. 

3.1  Abstract  formulation:  discussion 

We  consider  the  primal  'fine'  problem:  Find  U  £  U  such  that 

A(U;v)  =  F(v)  VueV,  (3.1) 

where  A  :  U  x  V  — >  M  is  a  semilinear  form  and  F  :  V  — >  R  is  a  linear  functional.  Typi¬ 
cally,  this  is  a  variational  form  of  a  highly  complex  physical  problem  involving  a  nonlinear 
operator.  A,  and  a  sought  solution,  U. 

Accepting  that  this  is  too  computationally  demanding  to  work  with  we  introduce  the 
following  'coarse'  physical  approximation  to  this  problem  as:  find  u  &Uq  C  U  such  that, 

A0(u-,v)  =  F0(v)  Vv€V0cV,  (3.2) 

where  Aq  :  Uq  x  Vo  -h ►  R  is  a  semilinear  form  and  Fq  :  Vq  —■ ►  M  is  a  linear  functional.  The 
resulting  modelling  error  is  then  U  —  u. 

For  practical  purposes  we  also  need  a  computational  formulation.  This  is:  find  Uh  £ 
K  C  UQ  such  that, 

A0(uh',v)  =  F0(v)  VveV0h.  (3.3) 

This  introduces  a  discretization  error:  u  —  U}lf  and  formed  the  subject  of  the  first  phase  of 
this  project. 

Let  J  :  U  — >  M  be  a  quantity  of  interest  (Qol)  then  the  total  error  in  the  computational 
model  is, 

J(U)  -  J{uh )  =  J(U)  -  J(u)  +  J(u)  -  J(uh)  .  (3.4) 

v -  v 

modelling  error  discretization  error 
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The  bulk  of  this  report  is  concerned  with  the  modelling  error  term  but  for  complete¬ 
ness  we  also  include  some  material  related  to  discretization  error  described  in  earlier  re¬ 
ports  as  some  of  the  description  is  independent  of  whether  or  not  the  error  is  due  to  mod¬ 
elling  error  or  due  to  discretization  error. 

3.2  Abstract  formulation:  an  outline  of  the  procedure 

In  subsequent  sections  we  describe  in  some  detail  various  representations  of  J(U)  —  J(u ) 
that  can  be  used  to  estimate  this  quantity  without  explicitly  computing  the  fine  solution 
U.  In  summary  an  outline  involves  the  following. 

1.  We  solve  (3.3)  to  generate  the  coarse  approximation  Uh  ~  u  «  U. 

2.  We  then  use  Uh  as  data  in  what  is  called  our  dual  problem  which  in  the  implementa¬ 
tion  described  here  involves  the  following:  Find  z  &Uq  such  that 

A'(uh]  v,  z)  =  J'(uh ;  v)  Vv  e  U0,  (3.5) 

where  A'(.  \ .)  is  the  Gateaux  derivative  of  .4(.; .)  (see  later  for  a  definition)  and 
where  Uq  is  an  appropriate  space. 

3.  We  then  compute  the  right  hand  side  in  the  estimate 

J{U)-J{u)^F(z)-A{uh-z).  (3.6) 

In  this  way  the  form  A{.\ .)  describing  the  fine  problem  is  needed  to  get  the  residual 
term  in  step  3  and  also  to  describe  the  linear  problem  in  step  2  but  note  that  the  only 
nonlinear  problem  that  is  solved  is  in  generating  Uh  in  step  1. 

The  above  procedure  partly  fits  the  framework  described  in  [3]  with  the  main  dif¬ 
ference  being  that  in  this  approach  we  use  A' ., .)  derived  from  ,4(.; .)  instead  of  some 
simpler  approximation  to  it  as  is  allowed  in  [3].  In  this  way  we  do  about  the  best  possible 
in  the  estimate  in  step  3  based  upon  only  having  u/,  available  with  the  cost  of  doing  this 
being  in  the  difficulty  of  setting-up  and  solving  (3.5). 

3.3  Abstract  formulation  detail:  Introductory  comments  about  the  error  J(U)  — 
J(u) 

In  this  section  we  give  some  expressions  to  approximate  the  difference  J(U)  —  J(u)  in 
terms  of  a  function  z  satisfying  a  derived  problem  (i.e.  a  dual  problem)  in  each  case.  The 
approximating  expression  that  we  use  is  actually  the  same  in  each  case  and  is  of  the  form 

J{U)  -  J (u)  «  F(z)  -  A(u ;  z)  (3.7) 

with  the  difference  between  the  cases  being  in  how  2  is  defined.  The  first  derivations  that 
we  consider  are  useful  when  the  approximate  solution  u  is  obtained  using  the  same  form 
A{.: .)  as  is  used  to  characterize  the  exact  solution  U,  i.e.  this  is  the  case  when  we  only  have 
discretization  error  to  consider.  The  problems  involving  2  involves  first  Gateaux  deriva¬ 
tives  of  A(.; .)  and  of  the  Qol  J(.).  For  the  derivation  which  is  needed  for  the  approxi¬ 
mation  which  we  use  when  we  also  have  modelling  error  we  additionally  need  higher 
Gateaux  derivatives  and  thus  these  are  described  before  the  derivation  is  done. 
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3.4  Representations  suitable  when  we  only  have  discretization  error 

Let  V  be  a  function  space,  let  J  :  V  — >  R  be  a  functional  and  let  U  £  V  and  u  £  V  be 
functions.  Also  let  eu  =  U  —  u.  The  representations  described  here  are  partly  derived  by 
just  using  the  fundamental  theorem  of  calculus  as  follows.  We  let 

g:[0, 1] M,  g(s)  =  J(u  +  seu)  (3.8) 


so  that 


J(U)  -  J{u)  =  g{  1)  -  g{ 0)  =  f1  g’(s )  ds.  (3.9) 

Jo 

This  leads  us  to  introduce  first  Gateaux  derivatives  of  functionals  which  are  defined  by 


v)  :=  iim  J("  +  z  -M  =  A.  J(U  +  sv] 

s^o  s  ds 


s=0 


so  that  with  this  notation  we  have 


J(U)  —  J(u)=  /  J'(u  +  seu;  eu)  ds. 
Jo 

We  similarly  define  a  Gateaux  derivative  of  A{.  \ .)  by 


A'(u;  v,  w)  =  —A(u  +  sv ;  w) 

ds 


s=0 


and  thus  we  similarly  have 

A(U\w)  —  A(u;w)  =  /  A'(u  +  seu;eu,w)  ds. 
Jo 


(3.10) 


(3.11) 


(3.12) 


(3.13) 


In  (3.11)  and  (3.13)  we  have  differentiation  in  the  direction  of  the  error  eu  which  in  the 
applications  is  unknown.  To  proceed  further  we  consider  differentiation  in  all  possible 
directions  v  and  seek  a  function  z  £  V  such  that 


/  A'(u  +  seu;v,  z)  ds  =  /  J'(u  +  seu;  v)  ds  for  all  v  £  V. 
Jo  Jo 


(3.14) 


If  we  can  obtain  such  a  function  z  then,  with  v  =  eu/  we  have  the  representation 

J(U)  -  J{u)  =  A(U ;  z)  -  A(u;  z)  =  F(z)  -  A{u-  z).  (3.15) 


The  representation  (3.15)  with  z  defined  in  (3.14)  cannot  be  used  in  the  computations 
as  given  because  the  equation  for  z  depends  on  the  exact  solution  U  as  the  integral  involves 
a  path  from  u  (when  s  =  0)  to  U  (when  s  =  1).  To  get  something  which  can  be  used  in 
computations  we  must  replace  (3.14)  by  something  which  only  involves  terms  which  can 
be  computed.  As  it  is  assumed  here  that  we  can  obtain  u,  the  simplest  approximation  is  to 
use  a  one  point  quadrature  on  the  integral  involving  evaluation  only  at  s  =  0  giving 

A'(u-,  v,  z)  =  J'{u ;  v)  for  all  v  £  V  (3.16) 

with  now  for  this  z  £  V 

J{U)  -  J(u )  «  F(z)  -  A(u ;  z).  (3.17) 
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In  situations  when  we  can  obtain  an  estimate  eu  «  eu  we  can  do  a  bit  better  by  approxi¬ 
mating  (3.14)  by  replacing  the  unknown  eu  by  e„  and  then  approximating  each  integral  by 
the  trapezoidal  rule  to  give  the  problem  of  finding  z  G  V/,  ,  where  l  /,  is  a  finite  dimensional 
space,  such  that 

A'(u  +  eu;v,z)  +  A'(u;v,z)  =  J'(u  +  eu;  v)  +  J'(u;  v)  for  all  vE  Vh-  (3.18) 
J(U)  —  J(u)  »  F{z)-A{v,z)  (3.19) 

In  practical  situations  when  there  is  only  discretization  error  to  consider  this  works  well 
provided  a  suitable  eu  &  eu  can  be  obtained  and  provided  the  function  z  obtained  from 
(3.18)  is  sufficiently  close  to  the  function  2  satisfying  (3.14).  Observe  that  if  u  G  I4  satisfies 
A(u ;  v)  =  F(v)  for  all  v  G  14  then  we  must  seek  2  from  a  larger  space  than  14  as  otherwise 
our  estimate  of  J(U)  —  J(u)  would  just  be  0.  We  comment  on  this  further  in  subsection  3.7. 


3.5  The  definition  of  Gateaux  derivatives  again  and  Taylor's  series 


As  indicated  above.  Gateaux  derivatives  are  concerned  with  differentiating  forms  and 
functionals  in  the  direction  of  functions.  Specifically,  with  V  denoting  a  function  space 
we  have,  repeating  the  definitions  above,  that 


J'(u;  v) 


A!{u\  v,  w) 


J(u  +  sv)  —  J(u) 

lim - 

s 


—~J(u  +  sv) 
as 


8= 0 


A(u  +  sv,  w)  —  A(u;  w) 

lim - 

s 


— A(u  +  sv,  w) 

as 


s= 0 


(3.20) 

(3.21) 


The  semi-colons  in  both  cases  indicate  that  the  terms  are  possibly  non-linear  in  the  terms 
to  the  left  of  the  semi-colon  and  are  linear  in  the  terms  to  the  right  of  the  semi-colon.  In  the 
case  of  A'(u;  v,  w)  the  derivative  is  respect  to  the  term  just  to  the  right  of  the  semi-colon  (i.e. 
v).  It  should  be  noted  here  that  in  general  if  we  swap  v  and  w  we  get  a  different  number 
although  when  quasi-static  problems  are  described  we  will  have  situations  in  which  we 
have  terms  which  are  symmetric  in  v  and  w. 

By  differentiating  these  first  derivatives  in  a  similar  way  we  get  second  Gateaux 
derivatives  and  we  can  continue  this  to  define  higher  derivatives,  i.e. 


J"(u;vi,v2) 

:=  4-J'(u  +  svi;v2) 

s=o 

5 

(3.22) 

J"'(u-V  l,v2,v3) 

:=  —J"(u  +  sv  r,v2,v3) 

as 

5 

s= 0 

(3.23) 

etc. 

(3.24) 

In  terms  of  these  higher  derivatives  we  can  write  Taylor's  series  as 


g(s)  =  J(u  +  sv) 


2  3 

5(0)  +  sg’{  0)  +  S-g"{  0)  +  S-g'"{  0)  +  •  •  •  (3.25) 

J(u)  +  J'(u ;  sv)  +  -  J"(u ;  sv,  sv)  +  -  J'" (u ;  sv,  sv,  sv)  +  •  -(3.26) 

Li  o 
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3.6  Preliminaries:  Taylor's  series  in  R  with  remainder,  the  Peano  Kernel  repre¬ 
sentation  and  the  Trapezoidal  rule 


Before  we  consider  another  representation  for  J(U)  —  J(u )  we  need  a  result  concerning 
the  error  in  the  trapezoidal  rule  which  we  can  obtain  as  follows. 

First  let  i,i€R  and  let 


(x-t)l  :  = 


(3.27) 


(x  —  t)n,  if  x  —  t  >  0, 

0,  otherwise. 

V  ' 

Then  to  obtain  Taylor's  series  with  the  remainder  in  integral  form  for  a  function  which  is 
ri  +  1-times  continuously  differentiable  we  consider,  for  0  <  x  <  1,  the  following  term 
which  we  re-express  using  integration  by  parts  repeatedly. 

1 


Rn  = 


™  Jo 


f\x  -  f)^/(n+1)(f)  d  t  =  -  Ax  -  t  )nf^+1\t)  d  t 

Jo  n-  Jo 


(3.28) 


1 


n! 


1 


=  Rn- 1  +  —  (x  -  t)nf(n)(t)  =  Rn- 1 - Txn/(n)( 0),  by  parts,  (3.29) 


J  o 


n\ 


n 

■■■  =  Rq  —  ^  — xfe/(fc)(  0),  by  induction. 


k= 1 


Now  since 


we  obtain  the  result 


Ro  =  /  /(i)  (t)  d t  =  f(x)  -  /( 0) 


f(x)  =  Y,^kfik\0)  +  Rr, 


k= 0 


(3.30) 

(3.31) 

(3.32) 


Next,  if  L  is  a  linear  functional  for  which  we  can  interchange  the  L{.)  operation  and 
the  integral  operation  over  [0, 1]  then  we  have 

n  i  1  pi 

L(f)  =  T,TiL(xk)f{k\°)  +  -  L((x-t)l)f(n+1\t)dt,  (3.33) 

k= 0  K'  U' 

where  in  the  integral  the  L{.)  operation  is  being  done  on  a  function  of  x.  This  representa¬ 
tion  is  particularly  useful  when  the  operator  L(.)  is  such  that  L(xk)  =  0  for  k  =  0, 1,  •  ■  •  ,  n 
as  then  the  representation  reduces  to 

L(f)  =  -  f1  dt.  (3.34) 

n'  Jo 

This  is  known  as  the  Peano  kernel  representation  of  L(f). 


We  now  apply  the  Peano  kernel  representation  to  the  following  operator  L  given  by 

L(f)  :=  I"  f(x)  dx  -  /(0)^/(1).  (3.35) 

Jo  z 

which  gives  the  error  in  the  trapezoidal  rule.  We  have  L(l)  =  I  Ax)  =  0  but  L(x2)  A  0  ar|d 
thus  in  the  Peano  Kernel  result  we  have  n  =  1.  With  0  <  t  <  1,  the  integrand  gives 

L((x  -  t)+)  =  (x  -  t)+  dx  - 

=  (*  -  t)  dx  -  =  \  (t1  ~  ~  I1  “  0)  =  ~  !)<- 
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3.7  Further  comments  about  the  order  of  magnitude  of  J(U)  —  J(u )  and  the 
accuracy  needed  in  the  approximation  of  z 

With  the  preliminaries  of  the  previous  subsections  it  is  now  convenient  to  return  briefly  to 
the  versions  of  z  considered  in  subsection  3.4.  Firstly,  suppose  that  z  is  defined  by  (3.14) 
giving  the  exact  representation  (3.15)  and  the  set-up  is  such  that  U  G  V  and  u  G  14  C  V  is 
of  the  form  u  =  Uh  =  U  —  eu,  (where  h  indicates  the  mesh  size)  and  satisfies 

A(uh',v)  =  F(v)  for  all  v  G  14.  (3.37) 

In  this  case  J(U)  —  J(uh)  is  zero  if  U  =  Uh  and  is  also  zero  if  z  G  V),.  These  two  different 
cases  when  J(U)  —  J(u )  is  0  enables  us  to  deduce  something  about  the  magnitude  of  this 
term.  Mathematically  then  for  any  Zh  G  V),  we  have 


J{U)  -  J(uh)  =  F(z)  —  A(uh;  z)  (3.38) 

=  F{z)  -  A(uh ;  z)  -  ( F(zh )  -  A(uh;  zh))  (3.39) 

=  F{z  -  zh)  -  A(uh ■  z  -  zh)  (3.40) 

=  F(z  -  zh)  -  A(U  -eu;z-  zh)  (3.41) 

=  F(z  -  zh)  -  A(U\z  -  zh)  +  A'(U;eu,z  -  zh)  4 -  (3.42) 

=  A’(U;eu,z-zh)  +  --  -  .  (3.43) 


With  A!{.  \ ., .)  being  linear  in  terms  after  the  semi-colon  this  tells  us  that  the  magnitude  of 
J ( U )  —  J ( Uh )  depends  on  the  product  of  the  magnitudes  of  eu  and  z  —  Zh ■  Thus  in  situations 
when  the  Qol  functional  </(.)  is  such  that  the  function  z  is  sufficiently  smooth  we  have 
typically  that  the  quantities  eu  and  z  —  Zh  are  0(hp)  in  magnitude,  where  p  depends  on 
the  degree  of  the  piecewise  polynomials  used,  and  J(U)  —  •/ (  «/, )  is  (D(h2p).  That  is,  we  are 
likely  to  approximate  the  Qol  more  accurately  than  we  approximate  U. 

The  function  z  defined  in  (3.14),  which  we  are  not  likely  to  be  able  to  compute,  is 
the  not  the  same  as  the  function  z  defined  in  (3.18)  which  involves  eu  instead  of  eu,  a 
finite  dimensional  space  V),  and  the  trapezoidal  rule  approximation  of  the  integral.  To 
distinguish  between  these  functions  let  z  €E  V,  Zh  €  V ,  Zh  €  V/,,  and  zjr  G  V),  satisfy 
respectively 


J'(u  +  seu\  v)  d.s 
J'(u  +  seu;  v)  d.s 
J'(u  +  seu\  v )  d.s 


for  all  v  G  V, 
for  all  v  e  V, 
for  all  v  G  V),, 


J\u  +  eu;v)  +  J\u ;  v )  for  all  v  G  V/,- 


(3.44) 

(3.45) 

(3.46) 

(3.47) 


The  function  2^  is  what  we  actually  compute. 
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We  next  compare  z,  Zh,  zy  and  Zh  and  comment  on  the  properties  that  we  need  from 
eu  and  14  to  produce  a  sufficiently  good  estimate  Zh  ~  z  for  the  purpose  of  the  error 
estimate.  In  the  reasoning  we  assume  throughout  that  u  and  U  are  not  too  close  to  critical 
points  such  as  limit  points  and  bifurcation  points  so  that  the  linear  problems  (3.44)-(3.47) 
have  unique  solutions  which  are  not  too  sensitive  to  changes  in  the  data.  We  will  also 
assume  that  ||eu||  is  0(hp).  For  the  estimate  of  the  error  we  will  assume  that  this  can  be 
done  by  constructing  eu  so  that  \\eu  —  eu\\  =  0(hv  11 ).  In  such  a  case  there  is  only  a  change 
of  0{Kp+l)  in  the  data  from  (3.44)  to  (3.45)  and  thus  || z  —  Zh\\  =  0(hp+1).  To  have  a  similar 
change  when  we  compare  (3.45)  and  (3.46)  requires  that  V),  is  sufficiently  large  that  this  is 
the  approximating  power  of  this  space.  That  is,  if  we  get  0(hp )  approximations  using  the 
space  14  then  we  need  14  D  14  to  contain  polynomials  of  one  degree  higher  to  have  an 
approximating  capability  of  0(hp+1).  When  this  is  the  case  we  have  \\zh  —  zy,  I  =  0(kp+1). 
Finally,  to  compare  (3.46)  and  (3.47)  we  note  that  by  using  (3.36)  we  have 

J'(u  +  seu;v)  ds  -  +  eu\v )  +  J'(u\v )) 

1 

J"'(u  +  seu ;  eu,  eu,  v)s(s  -  1)  ds 

u  +  seu- v,  z)  ds  -  ^{A!(u  +  eu;  v,  z)  +  A’{u\  v,  z)) 

1 

2  j  A"'(u  +  seu- eu,  eu,  v,  z)s(s  -  1)  d.s. 

With  ||en||  being  0[hp )  and  the  integrands  involving  products  of  such  terms,  the  trape¬ 
zoidal  rule  approximation  only  changes  the  data  by  0(h2p)  and  we  have  ||4,  —  Zh\\  = 
0{h2p).  Combining  all  these  comparisons  gives  us  that  \\z  —  zy  |  is  0(kp+1).  As  a  conse¬ 
quence,  when  we  compare 


J(U)  —  J(u)  =  F{z)  —  A(u ;  z)  with  F(zy)  —  A(u;  zy)  (3.50) 

we  have,  similar  to  the  derivation  (3.38)-(3.43),  that 

F(z)  -  A(u ;  z)  -  ( F{zh )  -  A(u;  zh))  (3.51) 

=  F(z  -  zh)  -  A(u ;  z  -  Zh)  (3.52) 

=  F(z  -  zh)  -  A(U  -eu;z-  zh)  (3.53) 

=  F(z-zh)  -  A(U-,z  -  zh)  +  A'(U;eu,z  -  zh) -] -  (3.54) 

=  A'(U;eu,z-  Zh)  4 -  (3.55) 


With  ||eu||  being  0(hp)  and  \\z  —  Zh\\  being  0(hp+1)  the  error  in  the  estimate  using  Zh  is 
0(h2p+1)  in  magnitude  which  gives  us  asymptotic  exactness  as  J(U)  —  J(u)  is  0(h2p). 

The  asymptotic  exactness  just  described  relied  on  eu  being  such  that  \\eu  —  eu\\  = 
0(kp+1)  and  that  V),  being  such  that  |  zh  —  5/,  |  =  0(hp+1).  We  can  relax  these  requirements 
at  the  expense  of  reducing  to  a  consistent  estimate  of  J(U)—J (  u)  if  it  is  possible  to  compute 
eu  such  that  \\eu  —  eu\\  <  ci||eu||,  where  c\  is  a  constant  satisfying  0  <  c\  <  1,  and  14  14 

is  sufficiently  larger  than  14  to  generate  a  sufficiently  better  approximation  than  would  be 
obtained  if  only  14  is  used. 


(3.48) 


(3.49) 
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3.8  Remarks  about  when  critical  points  occur 

The  above  assumes  that  we  are  not  close  to  critical  points  and  thus  to  complete  these 
comments  we  remark  as  to  how  such  points  are  characterised  in  this  abstract  setting.  Now 
in  what  has  been  given  so  far  we  have  that  U  £  V  satisfies 

A(U;v)  =  F(v )  for  all  v  £  V.  (3.56) 

In  the  application  described  in  this  report  this  brief  notation  hides  some  of  the  detail  as  in 
the  quasi-static  problem  to  be  considered  the  function  U  represents  a  displacement  field  at 
time  t  as  a  result  of  a  loading  P(t).  That  is  we  have  the  situation 

A(U(t),P(t)-,v)  =  F(v)  for  all  v  €  V.  (3.57) 

(In  our  application  F(.)  does  not  vary  with  t.)  Now  if  we  consider  an  increment  of  time 
At  and  let  AU  =  U(t  +  At)  —  U (t)  then  at  time  t  +  At  we  have 

F(y)  =  A(U(t  +  At),  P(t  +  At);  v) 

=  A(U(t),  P{t  +  At);  v)  +  A'(U(t),  P(t  +  At);  AU,  v)  + 

=  (A(U(t),P(ty,  v)  +  A  t¥LjtA(U(t),  P(t);v)  +  •  •  • ) 

+A'(U (t),  P(t);  AU,  v)  +  ■  ■  ■ 

By  using  (3.57),  dividing  by  At  and  letting  At  — >  0  gives 

BU  BP  B 

A\U(t),P(ty,-^,v)  =  -  ——A(U(t),P(ty,v)  for  all  v  e  V  (3.62) 

OTT 

which  is  an  equation  for  the  velocity  If  at  some  time  tc  the  form  A'(U (tc),  P(tc ); ., .)  is 
such  that  there  exists  a  function  w  €  V,  w  /  0  with 

A'(U(tc),  P(tc);  ip,  w)  =  0  for  all  ^  G  V  (3.63) 

then  we  have  what  is  known  as  a  critical  point.  In  such  a  case  we  do  not  have  a  unique 
solution  of  a  dual  problem  such  as  (3.16)  since  if  z  is  a  solution  then  z+/3w  is  also  a  solution 
for  all  /3  G  R.  When  we  are  at  such  a  point,  letting  a  =  tv  in  (3.62)  gives 

^0 tc)^A(U(tc),P(tcy,w )  =  0  (3.64) 

and  one  or  both  of  the  terms  in  the  product  must  be  0.  The  limit  point  case  corresponds  to 
when 

BP 

—  (tc)  =  0  (3.65) 

3  p 

and  thus  if  the  loading  has  been  such  that  >  0  for  0  <  t  <  tc,  i.e.  the  loading  has  been 
increasing,  then  it  becomes  necessary  to  decrease  the  loading  for  t  near  tc  with  t  >  tc  in 
order  to  follow  the  solution  path.  This  will  not  be  pursued  here  with  the  purpose  of  these 
comments  being  just  to  indicate  that  the  techniques  of  representing  the  error  J(U)  —  J(u) 
in  terms  of  a  dual  solution  z  breaks  down  when  critical  points  are  encountered. 


(3.58) 

(3.59) 

(3.60) 

(3.61) 
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3.9  A  general  expression  for  J(U)  —  J(u)  when  we  have  modelling  error 

We  now  consider  the  difference  J(U)—J(u)  again  without  specifying  how  the  approximate 
solution  u  is  obtained  or  how  the  function  z  is  obtained  and  as  such  it  is  suitable  in  cases 
where  we  may  solve  a  different  problem  to  generate  them,  i.e.  we  have  modelling  error. 
The  representation  given  here  is  similarly  derived  in  [3]  and  see  also  [2]  which  contains 
similar  results. 

As  before,  let  F(v)  be  a  linear  functional  defined  on  a  space  V  and  suppose  that  U 
satisfies  the  following 

A{U-  v )  =  F(v),  for  all  v  G  V.  (3.66) 

Related  to  (3.66)  we  also  suppose  that  Z  G  V  satisfies 

A'(U;  v,  Z )  =  J'(U ;  v),  for  all  v  €  V.  (3.67) 

As  before  u  «  U  and  we  let  z  m  Z  be  obtained  in  some  way  and  let  the  errors  in  these 
approximations  be  defined  by 

eu  =  U  —  u  and  ez  =  Z  —  z.  (3.68) 

Let 

go  (s) 

9i(s) 

and  observe  that 

go(l)  =  J(U),  go(0)  =  J(u)  +  F(z)-A(u-,z),  gi(l)  =  0.  (3.71) 

Thus 

J{U)  -  J{u)  =  <7o(l)  -  go(0)  +  F{z)  -  A{u-  z).  (3.72) 

Now  by  the  trapezoidal  rule  with  remainder  we  have 

5o(l)  —  <?o(0)  =  j  5o(s)ds  =  ^(s-o(O)  +5o(l))  +  s(s  -  l)g'o  (s) ds.  (3.73) 
For  the  derivative  of  go  we  have 

g'o(s )  =  J'(u  +  seu;eu)  +  F(ez)  -  A'(u  +  seu;eu,z  +  sez)  -  A{u  +  seu\ez)  (3.74) 


=  9i  0' 

)  +  F(ez 

)  -  A(u  +  seu;ez). 

(3.75) 

Thus 

m  = 

0,  g'o(0)  =  gi(0)  +  F(ez)  -  A(u;ez) 

(3.76) 

and  we  have 

3o(l)  -5o(0) 

1 

2 

(. g'o(° )  +  Jo  s(s~  1)ffo/(s) ds) 

(3.77) 

1 

2 

(0)  +  F(ez)  -  A(u;ez)  +  J  s(s-l)g'o(s)  dsV 

(3.78) 

=  J(u  +  seu)  +  F(z  +  sez)  -  A(u  +  seu\  z  +  sez),  (3.69) 

=  J'(u  + seu-eu)  -  A'(u  + seu-eu,z  + sez),  (3.70) 
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As  <71(1)  =  0  we  have 


9i  (0)  =  -flg[{s)ds  (3.79) 

Jo 

=  -  (J"(u  + seu\eu,eu)  -  A" {u  + seu\eu,eu,z  +  sez) 

Jo 

~  A'(u  +  seu ;  eu,ez))  ds  (3.80) 

=  ~  (J"(u  + seu-,eu,eu)  -  A"(u  +  seu;eu,eu,  z  +  sez))  ds 

Jo 

+  {A{U-ez)-A{u-ez)).  (3.81) 

Combining  (3.78)  and  (3.81)  and  noting  that  F(ez)  =  A(U\  ez)  gives 
5o(l)  -  ffo(0)  =  F(ez)  -  A(u ;  ez) 

A"(u  +  seu\eu,eu,z  +  seu)  -  J"{u  +  se„;  eu,  eu)  +  s(s  -  l)g'o(s)  ds^  . 

Finally  combining  with  (3.72)  gives 

J{U)  —  J(u )  =  ( F(z )  —  A(u]  z))  +  ( F(ez )  —  A(u ;  ez))  +  smaller  terms.  (3.82) 

Again  this  is  back  to  the  estimate 

J{U)  -  J(u)  w  F(z)  -  A(u;  z )  (3.83) 

and  the  simplest  approximation  z  that  we  can  obtain  assuming  that  only  u  is  available  is 
to  take  the  function  z  satisfying 


A'(u;  v,  z)  =  J'(u]  v )  for  all  v  G  V. 


(3.84) 


4  The  physical  problem:  Application  of  the  theory  to  membrane 
inflation 

In  this  section  we  consider  the  physical  models  in  which  we  are  going  to  apply  the  theory 
of  representing  the  error  in  a  Qol.  Specifically,  we  consider  the  equations  describing  the 
inflation  of  a  hyperelastic  membrane.  With  U_  denoting  the  exact  displacement  and  V 
denoting  the  exact  velocity  the  weak  forms  will  be  written  in  the  form 

)  =  F(  )  for  all  appropriate  ip  and  9  (4.1) 

in  the  dynamic  case  and  in  the  form 

a(LLi  VO  =  FW  f°r  aH  appropriate  ip,  (4.2) 

at  each  time  t,  in  the  quasi-static  case.  In  this  latter  case  the  form  a(U_\  ip)  is  sometimes 
also  written  as  a(U]  ip) a,  when  we  wish  to  explicitly  show  the  domain  involved,  and  can 
also  be  written  as  a(U_,  d)o,  when  we  also  wish  to  explicitly  indicate  that  the  form 
depends  on  a  loading  term  P(t)  at  time  t. 
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4.1  The  inflation  of  an  elastic  membrane 


We  assume  that  the  reference  region  of  the  body  being  deformed  is  given  by 

n3D  =  {(X1,X2,X3):  (X1,X2)en,  \X3\<h0/2}  (4.3) 

and  we  are  considering  the  deformation  over  a  time  interval 

0  <  t  <  T.  (4.4) 

The  components  X\,  X2  and  X3  are  with  respect  to  a  fixed  cartesian  coordinate  system.  We 
suppose  that  \ho\  is  small  as  part  of  a  membrane  idealization  of  how  the  sheet  deforms.  In 
the  axi-symmetric  case  FI  is  a  circle  and  if  the  circle  has  radius  1  then  we  can  take  ft  = 
(0, 1)  =  {0  <  r  <  1}  with  respect  to  cylindrical  polar  coordinates  r,  'tj>,  x3.  Here  ft  is  the 
region  of  the  undeformed  mid-surface. 

The  deformation  of  the  body  is  described  by  the  mapping 

X-+  x:=  X  +  U,  (4.5) 

where  x  =  x(X,  t )  is  the  deformed  position  and  u{X_,  t)  is  the  displacement.  The  deforma¬ 
tion  gradient  F3/3  (which  is  the  Jacobian  of  this  mapping)  is  given  by 

F3D  =  I3D  +  V3DU,  where  V3dC/=(|^.,|^,K(-).  (4,6) 

Here  I, 3 d  is  the  identity  tensor.  Note  that  in  terms  of  components  F373  and  V3 dLL  are 
represented  by  3  x  3  matrices.  For  our  membrane  model  we  use  the  notation  F  and  XU 
for  the  corresponding  3x2  quantities 


/  Fn  Fi2\ 

F  =  I  F-2i  F22  I 

\F:u  F32J 


and  XU 


(  dU  dU  \ 

\mcu  dx2)  ■ 


(4.7) 


We  restrict  consideration  to  incompressible  materials  which  means  that 


det  F:!/)  =  1,  and  the  density  p  =  const. 


(4.8) 


The  simplifications  to  this  when  we  have  an  axi-symmetric  problem  and  we  use  cylindrical 
polars  is  that  the  deformation  of  the  mid-surface  is  given  by 

/ 1  +  U[  0  \ 

(r,^,0)  ->  (r  +  Ui,ip,U3)  giving  F=  0  1  +  ^  (4.9) 

\  U’3  0  J 

with  '  denoting  differentiation  with  respect  to  r. 

Let  er  denote  the  Cauchy  stress  tensor  and  let 

LI  3D  =  (det  F3jD)F3^ct  =  F“^cr  (4.10) 

denote  the  nominal  stress  tensor.  Let  n  denote  the  unit  normal  to  the  deformed  mid¬ 
surface.  This  direction  corresponds  to  the  x:j-di  recti  on  in  the  reference  configuration.  In 
the  membrane  approximation  of  how  a  thin  sheet  deforms  it  is  assumed  that  an  =  0  which 
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corresponds  to  the  fact  that  the  stress  components  corresponding  to  the  direction  n  in  an 
actual  thin  sheet  are  generally  much  smaller  in  magnitude  as  compared  to  the  tangential 
components.  This  in  turn  implies  that  the  third  column  of  II \D  is  0  with  respect  to  cartesian 
coordinates  and  as  a  consequence  the  membrane  first  Piola  stress  can  be  taken  as 


/fin  n2A 

nr  =  n12  n22  (4.U) 

\iii3  n23/ 

and  we  have 

II3D  :  V3dV^  =  II  :  V0,  (4.12) 

where  here  :  denotes  the  double  dot  product  operation,  i.e.  A  :  B  =  rj  AVJ  BtJ  = 
tr(A7  B).  With  this  result  and  when  |/iq|  is  small  we  have 


II3Z)  :  V3d0  dX. 


ho 


:  V0  d,Y,dA'2. 


(4.13) 


We  use  this  in  the  weak  form  of  the  equations  of  motion. 

Let  V  =  U  denote  the  velocity  vector.  We  describe  the  equations  of  motion  in  the 
general  3D  case  in  partial  differential  equation  form  as  the  system 


pt 

V 


Div  Il3£),  where  Div  II 3^ 


U  . 


(d  nn 
1  TDC[ 
dU12 

.  dn13 


+ 

dn2i 

+ 

an31  \ 

+ 

dn22 

+ 

dn32 

,  (4.14) 

+ 

dn23 

+ 

dU33 

(4.15) 

To  complete  this  form  of  the  description  we  need  initial  conditions  U[X_,  0)  =  U°,  V_(X_,  0)  = 
V°,  boundary  conditions  (i.e.  on  Oil)  and  a  constitutive  model  relating  the  stress  II3 p  to 
F:i/).  For  the  computational  model  and  the  error  representation  formula  we  need  a  weak 
form  which  can  be  obtained  as  follows.  We  take  the  dot  product  of  (4.14)  with  0  and 
integrate  over  the  space  time  domain 


Q  =  n  x  (0,  T) 


(4.16) 


and  we  similarly  take  the  dot  product  of  (4.15)  with  6_  and  integrate  over  Q.  From  the 
terms  involving  '0  we  can  re-express  things  using  the  divergence  theorem  which  brings  in 
the  pressure  loading  term  P  =  P(t )  and  we  can  make  use  of  (4.13).  Assuming  that  0  =  0 
on  dO,  we  have 


Div  Il3£)  •  0  d.Y 


V  3D0  dA  +  P 


(/1  x/2)-0dAidA2, 


—ho 


V0dXidX2  +  P 


(/1  x/2)-0cLYidA2, 


where  /  and  /  are  the  first  two  columns  of  F  given  in  (4.7).  We  combine  the  terms  in¬ 
volving  0  and  the  terms  involving  0  into  a  single  expression  in  which  the  initial  conditions 
are  only  imposed  weakly.  In  the  following  we  collect  together  the  unknowns  U_,  V  and  the 
test  vectors  0,  6_  as 


(4.17) 
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In  the  expressions  that  follow  we  use  the  usual  inner  product  notation  (.,  ,)c  to  mean 


(f,g)a  ■=  [  fgdG, 

Jg 

where  in  our  case  G  =  fl  or  G  =  Q.  Our  weak  forms  are  as  follows. 

In  the  quasi-static  case  we  ignore  the  pV_  term  and  we  have 

a(U_'i  VOo  =  0,  for  all  appropriate  ip,  for  0  <  t  <  T 

where 

n  :=  aiGZ;  YOn  -  P(t)a2{Uj,ip)n 

and  where 

=  h{)  I  n7  :  S7ip  dXdi, 

Jn  ~ 

a2(U',±)  =  J  {f  ±  x  /2)  ■  ip  dXdt. 

In  the  full  dynamic  case  we  include  the  pV_  term  and  we  have 

M  (o)  ’  (f)  )  =  (f) )  ’  for  a11  appropriate  ip  and  9, 

where 


(4.18) 

(4.19) 

(4.20) 

(4.21) 

(4.22) 

(4.23) 


a{U]±)Q  +  ph0(V,  ip) Q 


+ph0  (l 

{U-V,0)q  +  m,  O),0)a  +  (Z(-,0 ),V0n)  , 

(4.24) 

O 

o 

-e 

,9)n  +  pho(V_°,ip)n 

(4.25) 

where  a(.,  ,)q  involves  the  same  expression  as  in  the  quasi-static  case  except  that  now 
integration  is  over  Q.  Specifically  we  have 

a(U\  ±)q  =  ai  (. U\  YOq  —  a2  (£;  ip)Q  (4.26) 


with 


ai  (IZ;V0q  =  [  f  a(U;  ip)  dXdt  and  a2(U;ip)Q  =  /  f  P(f)a2(t/;^>)  dX 
Jo  Jn  Jo  Jn 


dXdt  (4.27) 


/o  jo 

with  the  integrand  expressions  being 

T  T 


aiOZ;YO  =  ^on7  :  X±  and  d2(U;ip)  =  (/x  x  /2)  ■  ip_. 


(4.28) 


The  dependence  of  II  on  U_  for  a  hyperelastic  material  depends  on  the  form  of  the  strain 
energy  function  W  and  can  be  expressed  in  the  form 


it 


dW 

~d¥ 


(4.29) 
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with  an  appropriate  interpretation  of  the  meaning  of  the  partial  derivatives.  In  the  case  of 
an  axi-symmetric  deformation  and  an  isotropic  material  we  can  give  more  explicit  repre¬ 
sentations  of  the  a(U_]  V')  term.  In  this  case  we  have  W  =  W (Ai,  A2)  where 


Xi  —  (1  +  U±)~  +  U 3  and  A2  —  1  + 


t/i 


(4.30) 


are  the  principal  stretch  ratios.  Then  it  can  be  shown  that  the  principal  stresses  are  given 
by 

dW  .  aw 


al  ~  T— 

o\\ 

To  simplify  the  notation  a  little  we  let 


and  (T  2  =  A2 


d\-2 

d2W 


(4.31) 


dW  dW  d2W  d2W  d2W 

dAi  dX2  dXf  dX2  dXidX2 


(4.32) 


Using  this  notation  it  follows  that  the  nominal  stress  components  are  given  by 


Ifn  =  t^(1  +  u[)  =  -7— (1  +  u[),  n  13  =  Sc/3  =  ^3,  n22  =  ^  =  W2.  (4.33) 

At  Ai  At  Ai  A2 


Then  by  noting  that  dX  =  rdr  and  letting  U_  =  (U\  .  t/3) 7  and  V  =  (V\ ,  V>,)T  we  have 


ai(m±)  =  ho  0^((1  +  UM  +  U’M  +  Wj±)  ,  (4.34) 

d2(U^)  =  (l  +  ^(-^i  +  (l  +  ^3).  (4.35) 


To  complete  the  description  and  to  add  some  rigor  we  need  to  indicate  the  function 
spaces  for  U_,  V  and  ij),  r  given  in  (4.17)  and  (4.23)  and  similarly  we  need  to  give  the  spaces 
for  the  test  vectors  and  6.  Note  that  we  have  first  spatial  derivatives  of  the  displacement 
U_  but  not  for  the  velocity  V  and  we  have  first  time  derivatives  of  both  U  and  V_.  For  the 
test  vectors  we  have  first  spatial  derivatives  of  i/j  but  there  are  no  spatial  derivatives  of  9 
and  there  are  no  time  derivatives  on  either  function.  Hence,  as  in  [1,  p.264],  the  unknowns 
and  test  vectors  are  such  that 

VeH1({0,T),L2(n)),  (4.36) 

^GL2((0,T),i/1(U)),  0eL2((O,T),L2(fi)).  (4.37) 

The  functions  V’  vanish  on  the  part  of  Oil  where  Dirichlet  boundary  conditions  are  im¬ 
posed. 

As  a  final  comment  about  imposing  the  initial  conditions  weakly  we  note  that  the 
displacement  initial  conditions  that  we  consider  are  such  that  U_°  =  0  or  U_°  corresponds 
to  a  uniform  pre-stretch  are  in  fact  satisfied  exactly. 


4.2  oi(.; .)  and  o2(.; .)  as  derivatives  of  functionals 

The  term  d\{U_]  'if))  and  the  integral  over  (0,1)  of  a2(lk  y)  defined  in  (4.28)  are  in  fact 
Gateaux  derivatives  of  appropriate  functionals,  i.e.  they  are  of  the  form 

rai(U;  VO  =  g[(Uji  VO  and  [  d2(U;,  VO  rdr  =  f  <j2 (U;  V0r  d'r  (4.38) 
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for  certain  functionals  g\  and  g2  .  Asa  consequence,  when  we  consider  the  dual  problem 
and  when  methods  for  solving  the  dynamic  problem  are  considered  we  have  to  consider 
terms  of  the  form 


ra[  (U_;  a,  ip)  =  g'[ ( U_ ;  a,  ip),  and  f  a2  (U_'i  Y0  rdr  f  g2  (U_;  a,  ip)  rdr  (4.39) 
—  —  Jo  Jo  ~ 

which  are  thus  symmetric  in  a  and  ip  and  as  a  consequence  we  can  construct  our  algo¬ 
rithms  in  such  a  way  that  each  linear  system  that  we  have  to  solve  has  a  symmetric  matrix. 

The  functionals  g\  and  g2  for  the  axisymmetric  case  are  as  follows. 

gi(U)-=rh0W  and  g2{U)  :=  ^(r  +  Hi)(-(r  +  UM  +  (1  +  Ui)U3).  (4.40) 

We  can  verify  directly  that  the  first  Gateaux  derivatives  of  these  do  indeed  give  the  expres¬ 
sions  for  ai(.; .)  and  02(4  •)  as  follows. 

We  first  consider  the  g\  term  and  we  use  the  notation  here  that  U_(s)  =  U_  +  sip  where 
ip  is  the  function  which  we  differentiate  in  the  direction  of.  We  similarly  let  Ai(.s)  and  A2(s) 
be  the  stretch  ratios  that  we  get  with  U_(s).  We  need  to  differentiate  with  respect  to  s  and 
evaluate  the  derivatives  at  s  =  0.  This  gives 


a'iGZ;  YO 


dAi 

ds 


s=0 


(1  +  U'i)ip'i  +  U31P3 
Ai 


and  A2  (JJj,ip) 


dA2 

ds 


s=0 


(4.41) 


Then  as  W  =  W( Ai,  A2)  the  chain  rule  gives 


g'l  m  ±)  =  rh0  (W1  a; (U;  l)  +  W2  y2  (U;  ±))  (4.42) 


and  we  get  the  term  ra2{Uj,  ip),  see  (4.34). 

Next  we  consider  the  g2  term  which  is  connected  with  the  pressure  loading  and  again 
use  the  notation  that  U_(s)  =  U_  +  sip  where  ip  is  the  function  which  we  differentiate  in  the 
direction  of  so  that 


3^2 (U{s))  =  (r  +  Ui  +  sipi)(-(r  +  Ui  +  sipi)(U %  +  sip'3)  +  (1  +  Ux  +  sip[)(U3  +  sip3)).  (4.43) 
Differentiating  with  respect  to  s  and  evaluating  at  s  =  0  gives 

^2(mt)  =  M-(r+UM+(l+U[)U3)  +  (r+U1)(-ip1U!i-(r+U1)ip,3+iP'1U3+(l+U[)ip3). 

(4.44) 

It  is  the  integral  of  this  over  the  space  domain  that  we  need  to  consider  and  we  can  re¬ 
express  things  using  integration  by  parts  on  the  ip[  and  ip3  terms.  Observe  that  since  ipi 
and  ip3  are  zero  at  r  =  1  we  have 


and 

Hence 
3 


-  f  (r  +  Ui)2ip3dr  =  2  f  (r  +  Ui)(l  +  U[)ip3  dr 
Jo  Jo 

[\r  +  UjUatt  dr  =  -  [\ (r  +  UM  +  (1  +  U[)U3)iPt  dr. 
Jo  Jo 


(4.45) 


(4.46) 


f1  g2(U;±)  dr  =  3  f1  (r  +  Ui)(—U3ipi  +  (1  +  U[)ip3)  dr  =  3  ['  d2(U;±)  rdr.  (4.47) 
Jo  Jo  Jo 
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4.3  The  Gateaux  derivative  A'(.\ ., .)  for  the  dual  problem  obtained  from  the 
dynamic  problem 

For  the  dual  problem  in  this  context  we  need  the  Gateaux  derivative  of  the  following  in  a 
direction  such  as  \~],  where  a  =  [  1  Y  d  =  (Y  | .  With  this  notation  we  have,  in  vector 

\flj  V«3  J  ~  \fhj 

notation, 

A((y);(f))  :=  a&±h  +  Pho(v,±h 

+pho  ((U-V,e)Q  +  (tf(,O),0)n  +  (V(;0),±h)  ,  (4.48) 

which  gives 

A((^) +5(?)  ;  (^))  :=  a^—  +  s^’th  +  Pho(V  +  sP,±)Q 

+ph0  (fSL-V  +  s(a-p),0)Q  +  (U(.,  0)  +  sa(.,  0),  0)n 
+m,0)+sP(.,0),±h).  (4.49) 

Thus 


A' (U- a,  f)Q  +  ph0(P,±)Q 

+ph0  {(a-P,0)Q  +  (a(,  O),0)a  +  (P(.,0),±h)  (4-50) 

A'(U]a,^)Q  -  phoiP^Q 

+ph0  (~(a,0)Q  -  %Q.h  +  (a(,T),0)n  +  (P(.,T),f)n)  (4.51) 


In  the  first  version  in  (4.50)  there  are  time  derivatives  on  ~  whilst  in  the  second  version 


in  (4.50),  which  was  obtained  by  integration  by  parts  in  time,  there  are  no  time  derivatives 
on  .  This  is  the  form  of  A'(.; .)  that  we  need  in  the  dual  problem  for  finding  (/>  and  9 
such  that 


>,fU\fa\  ($ 


vj  ’  \p  ’  \e 


)  =  A 


.  fa 


v  ’  V/?. 


)  for  all  suitable 


(4.52) 


The  "for  all  suitable  a  and  d"  condition  in  particular  requires  that  the  components  of  a(r,  t ) 
and  P(r,  t)  vanish  wherever  the  components  of  U(r.  t )  or  V(r,  t)  have  a  Dirichlet  boundary 
condition  and  hence  in  this  application  we  have  ai(0,  t)  =  B\  (0.  t)  =  0  (corresponding  to 
the  axisymmetric  condition)  and  a(l,  t)  =  (3(1,  t)  =0  (corresponding  to  the  clamped  edge 
boundary  condition). 

A  general  form  of  a  suitable  J(.)  is  given  by 


’{expression  in  U_,  V_,  U',  Vj}  drdf 


(4.53) 


+  [  r {expression  in  U(.,T),V(.,T),  U'(.,  T),V'(.,T)}drdt.  (4.54) 

Jo 
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There  should  be  no  time  derivatives  in  the  integrand  expressions.  The  Gateaux  derivative 
of  such  a  J(.)  is  of  the  form 

=  j  J  (a  ■  J_a  +  /?  •  J_p  +  a'  ■  +  ft  ■  J_g,)  rd?’df  (4.55) 

+  [\a(r,  T)  •  J_a(r,  T)  +  0(r,  T)  ■  Jp(r,  T)  (4.56) 

Jo 

a'(r,  T)  •  JQ,  (r,  T)  +  £'(r,  T)  •  Jg,  (r,  T))  rdr.  (4.57) 

By  considering  the  a(.,T )  and  /?(.,  T)  terms  we  get  conditions  on  p  and  9  at  time  t  =  T 
which  are 

ph0(a(.,T),9(.,T))Q  =  [\a(r,T)-Ja(r,T)+a'(r,T)-Ja,(r,T))rdr,  (4.58) 

Jo 

ph0(P(.,T),l(.,T))n,  =  [  (0(r,T)  ■  J_g(r,  T)  +  0{r,T)  ■  T))  rdr  (4.59) 

Jo 

which  must  hold  for  all  suitable  a(r,  T )  and  (3(r ,  T).  By  considering  the  6(r\  t)  terms  for 
0  <  t  <  T  we  get 

—pho(}P_  +  9,  P)Q  =  I  f  p  ■  Jp  +  0_  ■  J_p,  rdrdt.  (4.60) 

Jo  Jo 

Similarly,  considering  the  a(r,  t )  terms  for  0  <  t  <  T  we  get 

o>'{UJiQL^Q  ~  pho(a,9)Q  =  /  /  a-Ja+g'-J_a,rdrdt.  (4.61) 

do  Jo 

The  equations  (4.58),  (4.59),  (4.60)  and  (4.61)  give  a  linear  problem  for  -ip  and  9  which  has 
a  structure  which  is  similar  to  the  problem  of  determining  the  displacement  U_  and  the 
velocity  V  although  with  some  key  differences.  For  the  dual  problem  we  have  final  time 
conditions  (i.e.  at  t  =  T)  instead  of  initial  conditions.  Equation  (4.60)  is  similar  to  the 
condition  U_  =  V  which  is  imposed  weakly.  Equation  (4.61)  involves  the  term  a'(.; .) 
which  is  encountered  in  numerical  schemes  for  approximately  solving  for  U_  and  V  by 
using  a  Newton  iteration.  Recall  from  the  discussion  in  section  4.2  that  a' (Up  a,  ip) q  is 
unchanged  if  we  swap  a  and  ip. 

4.4  Expressions  for  a\  (U:  a ,  ip)q  and  a'2(U_;  a,  i/j)q 

We  now  complete  the  details  of  the  dual  problem  by  giving  expressions  for  the  a[(.;.,.)Q 
and  (4 (.; .,  ,)q  terms  needed  in  (4.61). 

Recall  that  we  have  already  shown  that 

rd(U-f)  =  g[ (U;iP)  =  rh0  (W1  A [(U-f)  +  W2  A'  (U-f))  .  (4.62) 

Taking  now  the  Gateaux  derivative  in  the  direction  of  a  we  have 

ra[(lL<x,±)  =  rh0  {(WU\'M;  a)  +  W12X'2(U;  a))A/1(t/;V;)  +  WiA?(C7;a,^) 

+  (WMU-a)  +  FF22A2 (Up  a))\2(Uj  ip)  +  W2X'^(U-a,±)) 

=  rh0  {WnX'M-,  a)X[(U]±)  +  W22X,2(U;a)X,2(U;iP) 
+W12(X'1(U;a)X'2(U-,f)  +  X[(U;ip)X2(U-,  a)) 

+(WiXi(U;  a,  f)  +  W2X2(U_;  a,  f)))  .  (4.63) 
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(4.64) 


The  first  Gateaux  derivatives  of  the  stretch  ratios  are  given  in  (4.41),  i.e. 

W  (TT.  ,.\  _  (!  +  +  G3V4  W  (JTw  x  _  (1  +  U[)a[  +  C/gOg 

^1  \LLiy)  \  1  OO  \  5 

—  Ai  Ai 

A,2(G;V’)  =  — ,  A/2(G;a)  =  ^. 

—  r  r 


(4.65) 


For  the  second  Gateaux  derivatives  of  the  stretch  ratios  we  immediately  have  A" (G;  a,  V>)  = 
0.  For  the  second  derivative  of  Ai  we  have 


AiA'/o/ja,  vo + i + <*^3 


and  hence 


A'/O^VO  = 


a'lfti  +  “3^3  -  Ai  (C7; «) Ai  ( t/; 


Now,  from  (4.35)  we  have 


where 


a2(G;V()<3=  [  [  P(t)  d2(U;'f)  rdrdt, 

Jo  Jo 


ra2(U ;  V>)  =  (r  +  G1)(-[/'V’i  +  (1  +  ^)V»3) 


(4.66) 


(4.67) 


(4.68) 


(4.69) 


and  from  this  it  follows  that 


r®2  {UjiQL.iij)  =  {r  +  Ui)(-a 3^1  +  ot^  3)  +  a^-GgV’i  +  (1  +  U[)^3).  (4.70) 

It  is  the  integral  of  this  that  we  need  and  by  using  integration  by  parts  on  the  last  term  and 
noting  that  a\  and  V>3  vanish  at  r  =  1  we  have 


Flence 


[  aiip3(l  +  U[)  dr  =  -  f  («iV4  +  ap/hX?"  +  G,)  dr.  (4.71) 

jo  Jo 

■1  r1 

ra'2(U;a,^)dr  =  -  +  (r  +  17i) (0-3^1  +  V4«i))  dr-  (4-72) 

)  Jo 


Again  we  have  a  quantity  which  is  symmetric  in  a  and  V’  as  explained  previously. 


4.5  A  summary  of  the  dual  equation  in  the  dynamic  case 

To  summarize,  the  dual  problem  is  as  given  in  (4.52),  i.e.  it  is  the  linear  problem  of  finding 


—  1  such  that 


^  (I) ;  ft)  ■(!)>= j,<  (I) ;  ft) >  for  al1  sm,able  ft 


where  A'(  yyj  '■■■■)  can  be  written  as 

^((y)  ’  ft)  ’  ft)'*  =  -  Pho(^t)Q  -  ph0(P,±  +  o)Q 


+pho(a(.,T),6)n  +  ph0{/3(.,  T),  V>)o 


(4.73) 


(4.74) 
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and  «/'(.;.)  is  of  the  form 


io  Jo 


(ex  *  ■  /  ,  -)-  •  J_a  ~\~  ot_  •  •/,./  ~(-  (3  •  J at )  rd?’df 


+ 


[  (a(r,T)  •  JQ(r,r)  +  [3{r,T)  ■  Jg(r,T) 

Jo 

QL(r,  T)  ■  Ja,  (r,  T)  +  £'(r,  T)  •  J#  (r,  T))  rdr. 


(4.75) 


Note  that  A'(.; .)  does  not  actually  depend  on  the  velocity  V).  By  using  (4.63)  and  the 
subsequent  expressions  we  have 


a/(ZZ;£L  VOq  =  /  [  ras(Uj,a,ip)  drdt  (4.76) 

~  Jo  Jo  ~ 


with 

/M 

a:i(U]a,'ijj)  =  (a^oq,  V’i 

W 


and  with  the  symmetric  3x3  matrix  G  given  by 

G  =  /i0Gi  -  P(t)G-2 


where 


(  W22 

WV2(1  +  U[) 

Wl2  U-i, 

r2 

r  X[ 

r  XT 

G,  = 

Wn  (1  +  u[) 
r  Ai 

Wn(i±un+WM 

(wn 

Wi 

“AT, 

wl2u’z 
\  r  IT 

Wn%r  +  Wi  — 
Ai 

+  u[f 
A? 

and 


(4.77) 


(4.78) 


(4.79) 


(4.80) 


5  Implementation:  Computational  schemes 

In  this  section  we  describe  the  computational  schemes  used  for  the  quasi-static  problem, 
the  dynamic  problem  and  for  the  dual  problem  derived  from  the  dynamic  problem.  As 
already  indicated,  the  quasi-static  and  dynamic  problems  are  nonlinear  and  require  an  it¬ 
eration  and  a  Newton  iteration  is  used  here.  In  all  cases  we  have  a  spatial  discretization 
and  in  the  last  two  cases  we  have  a  space-time  discretization.  Time  levels  are  actually  also 
needed  when  we  consider  quasi-static  problems  in  possibly  two  different  ways.  They  pro¬ 
vide  the  means  for  generating  sufficiently  good  starting  points  for  the  Newton  iterations, 
i.e.  when  we  seek  the  solution  at  the  loading  level  P(tj)  at  time  t:l  we  can  use  the  solu¬ 
tion  at  the  previous  time  tj- 1  to  start  the  Newton  iteration  with  guaranteed  convergence  if 
P(tj )  —  P{tj- 1)  is  sufficiently  small  and  we  are  not  too  close  to  critical  points.  A  sufficient 
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number  of  time  levels  is  also  needed  in  cases  when  the  Qol  involves  the  velocity  and  the 
purpose  of  the  computation  is  to  predict  this  Qol  using  just  the  quasi-static  solution. 

As  we  have  explained  in  the  previous  sections,  the  procedure  of  constructing  an  ap¬ 
proximate  dual  problem  and  having  a  solution  Zh  to  such  a  problem  enables  us  to  be  able 
to  estimate  the  error  in  an  approximation  ,/(«/,)  to  a  Qol  J(U)  in  a  relation  of  the  form 

J(U)  -  J(uh )  «  F(zh)  -  A(uh;  zh).  (5.1) 

In  situations  when  we  only  have  discretization  error  we  can  make  further  use  of  z /,  to 
determine  how  we  should  adaptively  refine  the  mesh  in  order  to  create  an  approximation 
J{uh)  for  which  F(z)  —  A(uh:  Zh)  \  is  at  or  below  a  desired  tolerance.  We  describe  how  this 
is  done  in  the  case  of  the  quasi-static  inflation  problem  in  subsection  5.5. 

5.1  Notation  and  things  common  to  all  the  problems 

Much  of  the  notation  used  in  the  description  is  common  to  all  the  problems  and  thus  we 
give  first  the  notation  for  the  space  and  time  discretization  and  we  give  the  piecewise 
polynomial  form  for  the  vectors,  i.e.  u,  y_,  i/j  and  0,  in  the  computational  schemes  after  a 
discussion  about  the  spaces  that  the  exact  solutions  to  each  problem  are  in. 

5.1.1  The  mesh  points 

For  the  discretization  we  assume  that  there  are  N  +  1  time  levels 

0  =  to  <  h  <  t2  <  ■  ■  ■  <  tw  =  T  (5.2) 

and  a  fixed  mesh  of  M  space  elements  (V,;_i .  r,  ')  where 

0  =  r0  <  n  <  r2  <  ■  ■  ■  <  =  1.  (5.3) 

The  meshes  used  need  not  be  the  same  for  each  problem. 


5.1.2  Remarks  about  the  function  spaces  for  UjVjip  and  9 

In  section  4  we  derived  the  expressions  for  a(.; .),  -4(.; .)  and  A'(.; ., .)  without  too  much 
comment  as  to  the  spaces  in  which  the  various  terms  should  be  in.  In  fact,  for  the  nonlinear 
problems  you  cannot  specify  precisely  the  function  space  for  the  displacement  U_  until  you 
have  specified  the  form  of  the  strain  energy  function  W  =  W( Ai,  A2)  in  order  to  know 
which  powers  of 

Ai  =  ((l  +  t/()2  +  t/f  )1/2  and  A2  =  l  +  ^  (5.4) 

appear.  Also,  even  when  you  are  in  a  position  to  specify  a  function  space,  so  that  a(U_] .) 
and  A(U_] .)  are  properly  defined  for  all  U_  in  the  space,  there  is  still  likely  to  be  a  far  bit 
of  additional  analysis  needed  to  determine  if  such  a  problem  has  a  unique  solution  from 
the  space.  A  rigorous  mathematical  analysis  of  which  space  U\  and  U3  are  in,  to  the  extent 
of  determining  the  values  of  q  >  1  such  that  U-  €  Lq(Q),  does  not  actually  affect  how 
we  choose  the  finite  dimensional  spaces  used  in  the  computations.  It  is  just  the  order  of 


-25- 


the  space  and  the  time  derivatives  appearing  in  the  expressions  which  we  need  to  take 
account  of  in  addition  to  requiring  conditions  such  as  U\  (0,  t)  =  0  for  our  axisymmetric 
deformation. 

For  the  a(.; .)  term  appearing  in  (4.19)-(4.22)  we  require  that  the  approximation  u  ~  U_, 
with  u  =  u(r,  t),  satisfies 

G  Ff1(f2),  i  =  1,3  and  tti(0, f)  =  0  (5.5) 

and  also  satisfies  the  boundary  conditions  at  r  =  1.  By  considering  the  term  Vh  we  simi¬ 
larly  take  test  vectors  such  that 

ipi(-,t)  G  iF1(n),  i  =  l,3  and  ip^O,  t)  =  ipi(l,  t)  =  ^3(1,  t)  =  0.  (5.6) 

The  formulation  of  the  dynamic  problem  appearing  in  (4.23)-(4.25)  involves  space 
derivatives  and  time  derivatives  and  we  now  have  u  ~  C/  and  »  k  V.  We  take  these 
approximations  such  that 

Ui  G  Ff1((0,  T),  fl1(ff)),  €  Ff1((0,  T),  L2(fl)),  2  =  1,3  and  «i(0,  t)  =  2>i(0,  t)  =  0 

(5.7) 

and  also  that  they  satisfy  appropriate  boundary  conditions  at  r  =  1.  (The  notation  used  in 

(5.7)  is  standard,  meaning  that 

G  H1^),  €  L2(fi),  ui(r,.)€H1(  0,T),  vi(r,.)eH1(  0,T),  2  =  1,  3.) 

(5.8) 

For  the  test  vector  ^  in  this  context  we  essentially  have  the  same  situation  as  in  the  quasi¬ 
static  case,  except  that  we  now  have  to  also  consider  the  time  dependence,  and  we  take 

ipi  G  L2((0,T),H1(n)),  2  =  1,3,  and  ^i(0,  t)  =  -0i(l,  t)  =  ip3(l,  t)  =  0  (5.9) 

as  there  is  no  time  dependence  on  ip  in  the  A(.; .)  expression  in  (4.24).  The  test  vector  6  in 
(4.24)  is  just  concerned  with  weakly  imposing  the  condition  that  the  velocity  is  the  time 
derivative  of  the  displacement  and  as  the  expression  involves  no  space  or  time  derivatives 
we  just  need 

0i€L2((O,T),L2(fi)),  2  =  1,3  with0!(O,t)  =0i(M)  =03(M)  =  0.  (5.10) 

We  actually  restrict  to  functions  6t  G  /^((O-  T).  II 1  (Q))  c  L2((0, T),  L2(fi))  in  the  imple¬ 

mentation. 

The  test  vectors  ip  and  6  in  the  dynamic  problem  are  the  unknowns  in  the  dual  prob¬ 
lem,  which  is  derived  from  the  dynamic  problem,  but  they  appear  in  a  different  way  as 
the  expressions  were  obtained  using  an  integration  by  parts  in  time  so  that  the  expression 
that  we  use  for  A'(.; ., .)  involves  first  time  derivatives  on  ip  and  0.  Thus  when  we  solve 
for  ip  and  0  in  the  numerical  scheme  we  need  to  be  able  to  integrate  the  time  derivatives  of 
these  and  hence  we  take 

ipi<EHl(( OjT),#1^)),  2  =  1,3,  and  ipi(0,t)  =  ipi(l,  t)  =  ^3(1,  *)  =  0  (5.11) 

and  we  similarly  we  take 

0i  G  1T1((0, T),  W1  (fl)),  2  =  1,3,  and  0^0, t)  =  6i(l, t)  =  03(1, t)  =  0.  (5.12) 
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In  the  dual  problem  it  is  the  vectors  labelled  as  a  and  /3  which  have  the  role  of  the  test 
vectors  and  recall  that  in  the  derivation  these  are  concerned  with  the  directions  in  which 
we  change  U_  and  V  respectively  when  the  Gateaux  derivatives  are  obtained.  Observe  that 
there  are  no  time  derivatives  on  a  and  (3  in  (4.51).  It  thus  follows  that  we  can  take 

oti  €  ^2((0,  T),  i  —  1,3  and  ai(0,  t)  =  ai(l,  t)  =  0:3(1,  t)  =  0  (5.13) 

and 

/3j€L2((0,T),L2(fi)),  *  =  1,3  and  /?i(0,  t)  =  /?i(l,  t)  =  /%(1,  t)  =  0  (5.14) 

although  in  the  computations  we  actually  take  f3t  €  ^((0,  T).  H 1  (il))  c  -^2 ((0,  T).  /^(fl)). 


5.1.3  The  piecewise  polynomial  approximations 

We  take  piecewise  polynomials  for  the  approximations  u  ss  U_  and  v  ~  V  and  for  the 
approximate  solution  ip,  0  in  the  dual  problem.  We  have  first  time  derivatives  on  the 
vectors  u,  v,  V7  and  0  and  thus  we  take  functions  which  are  continuous  and  piecewise 
linear  in  t  with  respect  to  the  time  levels  0  =  to  <  ti  <■■■  <  tN-  We  similarly  have  first 
derivatives  with  respect  to  r  on  both  components  of  the  first  three  vectors  and  thus  for 
these  vectors  we  take  functions  which  are  continuous  and  piecewise  polynomial  in  r  with 
respect  to  the  mesh  points  0  =  ro  <  r\  <  ■  ■  ■  <  rAy  =  1.  We  also  take  the  same  form  for  0. 
For  some  detail,  we  let  Hi,  ■  ■  •  ,  Hm  denote  the  spatial  basis  functions  so  that  at  any  time 
level  t  =  tj  we  have 


M 


M 


u(r,tj)  =  uj(r)  :=  y(r,tj)  =  yj (r)  :=  ^y^H^r), 


(5.15) 


i=  1 


i=l 


and  similarly  when  we  consider  the  dual  problem 


M 


M 


?Kr’tj)  =  r ) :=  ESlW.  -Mr)  :='£S«i(r). 


(5.16) 


i= 1 


i=l 


(The  number  of  basis  functions  M  depends  on  the  number  of  elements  M  and  the  degree 
of  the  polynomials  used,  e.g.  M  =  M  +  1  for  linears,  M  =  2M  +  1  for  quadratics  etc..) 
Between  time  levels,  i.e.  t,_i  <  t  <  tj  we  assume  that 


n(r,  t ) 


y(r,  t ) 


V>(r,  t ) 
9(r,  t ) 


tj  ~  tj- 1 
tj~t 
tj  ~  tj- 1 
tj~t 
tj  ~  tj- 1 
tj  ~  t 
tj  ~  tj- 1 


-  V)  +  (J.  _tjt  . 

vi~l(r)  +  f  ^  vJ(r), 

+  (j7ZTj~)  ’ 

,^7^)  ^(r)- 


(5.17) 

(5.18) 

(5.19) 

(5.20) 
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The  time  derivatives  of  u,  v  and  L  are  defined  for  tj- 1  <  t  <  t  j  and  are  given  by 


u(r,t.)  =  - - - (vP  (r)  —  uJ  (5.21) 

tj  tj- 1 

v(r,t)  =  —— i- - (v? (r)  —  v7-1(r))  ,  (5.22) 

tj  tj- 1 

=  — — - - (V^(r)  -  ±J~1(r))  ,  (5.23) 

tj  tj- 1 

0(r,t)  =  — - - (0*'(r)  -0j-1(r))  .  (5.24) 

tj  ~  tj- 1 

These  time  derivatives  are  piecewise  constant  in  time. 


5.2  The  quasi-static  problem 

In  the  computational  scheme  for  the  quasi-static  problem  we  determine  the  displacement 
u3  (r)  at  each  time  level  which  satisfies 

a(]jp;ip)  =  F{i>)  f°r  all  V’  ^  *S 'h,  (5.25) 

where  5),.  is  an  appropriate  subspace  of 


The  only  test  vectors  which  are  not  used  are  those  attached  to  displacement  degrees  of 
freedom  which  are  known  corresponding  to  the  points  r  =  0  and  r  =  1.  (5.25)  gives 
a  nonlinear  system  in  the  nodal  parameters  of  yP  which  are  unknown.  Let  (ur)(k> ,  k  = 
0, 1,  2,  •  •  •  denote  the  functions  corresponding  to  each  Newton  iteration  in  the  attempt  to 
solve  (5.25)  and  we  write  instead  a((yP)tk\  P(tj)-,ip)  for  a((uJ")^;^)  to  make  clearer  that 
the  form  changes  with  j.  As  already  indicated  a  possible  starting  iterate  is  (W)1'0'  =  yP~l 
and  subsequent  iterates  are  determined  so  that 

a{(yp)[k\  P(tj);^)  +  a\(yp)(k\  P{tj);  (ui)(k+V  -  =  F{f)  for  all  Sh.  (5.27) 

This  is  a  linear  system  for  the  change  (yP)(k+1^  —  (wr)(k>  in  the  iterates  and  the  linear  system 
has  a  matrix,  the  Jacobian  matrix,  which  is  symmetric  because  of  the  symmetry  property  of 
a' {{vP)(k\  P(tj)\ ., .)  discussed  in  subsection  4.2.  If  this  Newton  iteration  fails  to  converge 
with  the  loading  P(tj)  then  we  can  change  tj  to  be  closer  to  tj- 1  at  which  we  already  have  a 
solution.  This  modification  of  reducing  the  time  step  will  eventually  lead  to  a  convergent 
iteration  provided  we  are  not  at  or  very  close  to  critical  points  as  at  critical  points  the 
Jacobian  matrix  is  singular. 


5.3  The  dynamic  problem 

The  approximate  dynamic  solution  at  time  tj  involves  determining  uP  and  y3  satisfying  an 
equation  of  the  form 

A((^),Pfe);(|))  =  F((|)),  for  all  Sf,  deSdh,  (5.28) 
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where  S is  an  appropriate  space  constructed  using  the  functions  Hi,  ■  ■  ■  ,  Hm  to  describe 
the  spatial  dependence.  As  the  test  functions  only  need  to  be  such  that  ip(r, .)  £  /-2  ( (0 .  T)) 
and  9(r, .)  £  L2((0,  T))  we  can  choose  them  to  be  non-zero  only  on  one  time  interval  at 
each  stage  and  in  this  way  the  equations  we  solve  at  time  tj  only  involve  the  solution 
already  obtained  at  one  previous  time  level,  i.e.  at  tj- 1. 

To  get  started  we  need  the  solution  at  the  initial  time  to  =  0  and  this  is  obtained  by 
projecting  the  initial  displacement  U_°  and  the  initial  velocity  V°  into  the  approximating 
space  in  the  sense  that  u°  and  v()  satisfy 

[  u?(r)-9_rdr  =  f  t/°(r)-#rdr,  /  y°(r)-iprdr  =  j  V^(r)-iprdr  for  all  ip, 9  £  S%. 

Jo  Jo  Jo  —  Jo 

(5.29) 

This  gives  linear  systems  involving  mass  matrices.  Note  that  for  a  simple  initial  condition 
such  as  a  pre-stretch  state  corresponding  to  U_°  =  (jr,  0)7  we  have  U_°  £  Sh  and  hence  in 
this  case  we  simply  have  u°  =  U_°.  In  some  of  the  examples,  this  projection  is  however 
necessary  to  generate  the  initial  velocity  vector. 

We  consider  next  how  to  determine  the  solution  at  time  tj  for  3  >  1-  Now  whatever 
choice  is  made  for  ip  and  0  we  have  equations  of  the  form 


pho{u  —  v)  ■  0_ rdrdt  =  0  for  all  appropriate  9,  (5.30) 

r  (ai(u(r,t)-,ip(r,t))  -  P(f)a2(n(r,f);^(r,  t))  +  ph0v(r,  t )  •  pJ)  drdt  =  0 
for  all  appropriate  ip.  (5.31) 


We  get  different  schemes  for  different  choices  of  the  time  dependence  of  9  and  ip  with  in 
all  cases  a  quadrature  needed  to  approximate  the  time  integral  in  (5.31).  We  just  describe 
one  possibility  here  in  which  we  take  V’  and  9_  to  be  constant  in  time  on  (t;_  i ,  tj)  and  for 
the  quadrature  in  time  we  use  the  mid-point  rule.  This  gives  us 


{tj-tj- 1)  /  [aipuJ  1/2',ip)  -  P{tj_i/2)a2{vP  i/x;^)  + 


1—1/2., 


ph  o 


tj  tj- 1 


(yi  —  yj  1 ) 


where 


h- 1/2 


^'-i  +  tj 
2 


•  ip)  rdr  =  0 

(5.32) 

(5.33) 


and  where 

i  i  M 

yj  1/2  =  - {vj^  +  vj)  =  -  1  +  yJk)Nk{r).  (5.34) 

To  be  able  to  obtain  an  equation  involving  only  y?  we  need  to  express  yj  in  terms  of  yP  by 
making  use  of  (5.30). 


As  already  stated,  our  choice  is  that  9  is  independent  of  t.  and  we  have 


J  f(yp  —  yj  7) 


Now  since 


-(tj  —  tj-i)(yj  1  +  yj)  j  ■  9_rdr  =  0  for  all  appropriate  0  £  5^. 

’  (5.35) 

(i i}  +  y )  -  :<i  (5.36) 
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it  follows  that 


1 


giving 


id  —  y?  1  =  -(tj  —  tj- i)(V  1  +  id) 
-( u  —  y d_1)  —  y 3-1 


v J  = 


tj  tj- 1 


and  the  acceleration  term  is  given  by 


7  7  —  1 

ir  —  7r 


t  j  —  tj  —  X  t  j  tj  —  l  y  tj  tj  —  X 


id  —  id  1 


—  v 


,1-1 


Substituting  into  (5.32)  gives,  after  dividing  by  tj  —  tj-\, 

ai(id_1/2;  V>)  -  P^.i/a)^^-1/2;^) 
2pho  (  up  —  yp~l 


+ 


3  3 


^ - yi-1)  -i/Ardr  =  0 

tj  tj  —  i  y  tj  tj— 1  /  J 


for  all  appropriate  ip  £  S%. 


(5.37) 

(5.38) 


(5.39) 


This  is  similar  in  form  to  (5.25)  in  the  quasi-static  case  except  that  here  the  pressure  is 
evaluated  at  the  mid-time  tj_ i/2  and  we  have  the  extra  inertia  term.  The  Newton  iteration 
that  we  use  to  solve  this  equation  is  also  similar  to  the  quasi-static  case.  If  we  write  the 
equations  in  the  form 

g{;up\  ip)  =  0  for  all  appropriate  ip  G  (5.40) 

then  the  Newton  iteration  for  solving  this  involves  iterates  computed  from 

g({yP){k);±)  +  g'({uj){k)-,  (yp){k+1]  -  (yj){k\±)  =  o.  (5.41) 


If  we  let  y  =  (id)(fc+1)  —  (y?)^  and  let  (id  l!2pk^  =  (up  1  +  (yP)^)  then  the  Gateaux 
derivative  of  g  is  given  by 

d((id)(fc);y,d)  =  jf r  (a'i{{uJ~1/2)(k)-,y,±)  -  P(ij_i/2)o2((«J_1/2)(fc);y,^)) 

+77~~~T2  H '  rdr  (5-42) 


which  should  be  compared  to  the  corresponding  term  a'((yp)^;y,ip)  in  the  quasi-static 
scheme.  The  procedure  thus  involves  solving  this  nonlinear  equation  for  id  by  Newton's 
method  and  setting 


As  in  the  quasi-static  case, 
metric  matrix. 


3- 1 


-(id  -  y?  1)  —  yp  1. 


(5.43) 


the  linear  system  associated  with  these  equations  has  a  sym- 
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5.4  The  dual  problem 


To  approximately  solve  the  dual  problem  (4.74)-(4.75)  we  must  already  have  an  approxi¬ 
mation  u  «  U_  available  which  in  this  application  could  be  an  approximation  to  the  solu¬ 
tion  of  the  dynamic  problem,  from  which  the  dual  problem  is  derived,  or  it  could  even  be 
an  approximation  obtained  by  neglecting  inertia  effects  throughout.  Note  that  as  the  the 
expressions  explicitly  involve  the  final  time  t n  =  T  the  solution  procedure  must  be  back¬ 
ward  in  time  involving  first  determining  ?/>(.,  T)  and  (]_(..  T )  and  then  determining  tj) 
and  for  j  =  N  —  1,  N  —  2,  •  •  •  ,1)0-  In  the  description  we  let  S®  denote  the  appro¬ 

priate  subspace  of 


span 


H, 


0 

Hi 


H2\  (0 

0  J  ’  \lhj  ’ 


Hm\  (  0 

0  J  ’  Um 


(5.44) 


in  which  each  function  ifd  (r)  and  03  (r)  lie. 


5.4.1  The  final  time  t  =  tu  =  T:  determining  ipN  (r)  and  (r 


Our  discretization  of  (4.74)-(4.75)  is  to  take 

(a,  0N)n  =  f  (a(r)  •  Ja(r,  T)  +  a'(r)  ■  J_a, (r,  T))  rdr,  for  all  a£Sh,  (5.45) 

Jo 

{P,±Nh,  =  [  (P-J0(r,T)+p'-j0,(r,T))rdr,  for  all .  (5.46) 

Jo 

Both  of  these  involve  linear  systems  with  mass  matrices.  This  step  is  similar  to  the  step  of 
obtaining  uP  and  c°  in  the  dynamic  problem. 


5.4.2  At  time  t  =  tj,  j  <  N:  determining  ?/4  (r)  and  (P(r) 


As  with  the  dynamic  problem  we  can  choose  the  test  vectors  ( a  and  (3  in  this  case)  to  be 
non-zero  only  on  one  time  interval  and  thus  we  have 


~ph  o 

and 


fb+i  r1 


o  \tj+ 1  tj 


rb+i  r1 


(V’J+1  —  tpi)  +  6  )  •  f3  rdrdf  = 


rb+i  r1 


(P  '  J./3  +  fJ'  '  J-ff )  rdrdt 
(5.47) 


d[  (u(r,  t)\a,  i/j)  —  P(t)a'i(u(r,  f);  a,  i/j)  —  ph^a  ■  8 )  rdrdt 
tj+ 1  r1 


It,  Jo 


{a-J^y  +  a-  J_ai)  rdrdt. 


(5.48) 


These  resemble  in  form  equations  (5.30)  and  (5.31)  in  the  dynamic  problem  and  we  again 
have  the  situation  in  which  different  choices  of  the  time  dependence  of  the  test  functions 
gives  us  different  schemes.  We  follow  a  similar  approach  to  that  used  in  the  dynamic 
scheme  and  we  present  one  scheme  here  corresponding  to  the  case  that  a  and  f3  are  both 
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constant  on  each  time  interval  (tj,tj+ 1)  and  we  approximate  the  time  integrals  by  quadra¬ 
ture  and  using  the  mid-point  rule.  This  gives 

(tj+i  ~  tj)  J  {a[(uJ+1/2;a, ip3+1/2)  -  P{tj+1/2)d'2(uJ+1/2- rdr 


>o 


~ph o  /  a  ■  ( 
-1 
'o 


1  (6j+1-0j 


rdr 


V  tj+ 1  tj 

=  (tj+i  —  tj)  /  rdr  for  all  a  G  5^*,  (5.49) 


where 


fj+l/2  - 


tj  +tj+i 


and  uj+1/2(r)  =  u(r ,  fj+i/2) 


and  where 


±3+1/2  =  +  ±3+1)  =  ^  ((^)fc  +  (V^1)/,-)  JVfc(r)  G  Sf . 


(5.50) 


(5.51) 


In  this  case  we  need  to  use  (5.47)  to  express  t V  in  terms  of  ip3  so  that  we  can  get  an  equation 
which  only  involves  ip3  and  the  scheme  that  we  get  depends  on  the  choice  of  /3  and  we  just 
restrict  to  the  case  that  /3  is  constant  in  time  on  (tj,tj+i).  Then 

—pho  J  (V+1  —  ip3)  +  fy+19 — —(03  +  03+1)^J  ■  fdrdr 

=  (tj+i~tj)  f  (/?  •  J_b{u3+1/2)  +  (3'  ■  J_3'(u3+1/2))  rdr  for  all  (3  G  Sk  .  (5.52) 
Jo  ~  ~  ~ 

To  get  the  connection  between  ip3  and  6_3  we  let 

f.  _f.  .  M-l 

y  =  (y+1  -  y )  +  J±L — !(&>  +  e3+1)  =  ykNk{r)  g  sf  (5.53) 

k=\ 

so  that 


— pho  f  ■  ardr  =  (tj+\  —  tj)  f  (/3  ■  J_p(u3+1/2)  +  (3'  ■  J_pi(u3+1/2))  rdr  for  all  (3  e  Sk  . 

Jo  ~  Jo  ~  ~  ~ 

(5.54) 

This  is  a  linear  problem  involving  a  mass  matrix  for  the  components  yl  which  are  un¬ 
known.  Once  these  are  obtained  we  have  the  functions  and 


Q3  =  _ej+ 1  + 


+ 1  t 


-  (yt  —  (y+1  —  y )) 


which  gives 


1 


tj- (-1  tj 

2 


(0i+i  _^) 


_  _L_(y  _  (y+ 1  -  y  ))^ 


tj+i  tj  V  tj+i  tj 

-2  ,  2 

- 2  y  +  - 

(tj+ 1 —  tj)  —  tj+i  —  tj 


9j+1  + 


(tj+i-tj)VX 


(y+1  —  ip3) 


(5.55) 

(5.56) 

(5.57) 

(5.58) 
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Substituting  into  (5.49)  gives  the  system  for  w3 ,  i.e.  we  have 

(tj+ 1  -  tj)  J  (fi,1(«J+1/2;a,^J+1/2)  -  P(tj+1/2)a!2(ul+1/2\a,jp+1/2))  rdr 

~2pK 1,  a '  ( jzFF)1’  +  (e~+'  +  rdr 

=  (tj+i—tj)  J  (a  ■  J_a(y?+1/2)  +  a'  •  J.a'(yJ+1/2)^J  rdr  for  all  a  €  S® .  (5.59) 

Collecting  together  the  'ip3  terms  on  the  left  hand  side  and  the  known  terms  on  the  right 
hand  side  gives 


^J+19  ^  J  (a[(yP+1/2;a,^3)  -  P(tj+1/2)d2(]F+1/2]a,pA)^  rdr 

+ - — — -  /  a -  ip3  rdr 

(tj+ 1  ~  tj)  Jo 

=  ~^+19  Jq  (ai(MJ+1/2;«,^+1) -JP(^+i/2)a/2(ff7+1/2;a,Vc7+1)) 


2p/i0 


»i 


fi 


/o 


a  •  y-7  rdr  +  2pho  f  a  ■  (  8_3+1  — 
»i 

'o 


7 - rV’",+1  I  rdr 

(fjr'  +  l  ~tj)~  J 


(tj- 1-1  fj)  do 

+(ij+i  -  tj)  I  •  J_a(uj+1/2)  +  a'  •  J_a'(uj+1/ 2))  rdr,  (5.60) 


for  all  a  £  S' [3 .  As  the  expression  on  the  left  of  the  equal  sign  is  symmetric  in  a  and  ip3 
we  get  a  linear  system  with  a  symmetric  matrix  to  determine  the  coefficients  of  'ip3 .  Once 
these  are  obtained  we  obtain  the  coefficients  of  93  using  (5.55). 


5.5  Using  the  dual  solution  to  adaptively  refine  a  mesh 

As  already  stated  a  few  times,  the  theory  described  in  this  report  leading  to  the  estimate 

J(U)  -  J(uh)  «  F(z)  -  A(uh ;  zh)  (5.61) 

of  the  error  in  approximating  a  Qol  J(U)  by  ./(u/t)  is  applicable  whatever  the  source  of 
the  error,  i.e.  discretization  error,  modelling  error  or  both.  When  there  is  modelling  error 
repeatedly  refining  the  meshes  used  to  obtain  Uh  and  Z],  becomes  pointless  as  F(zh)  — 
A(uh\Zh)  ■/>  0  but  instead  converges  to  the  modelling  error  alone.  In  this  situation  we 
would  just  wish  to  refine  until  two  approximations  of  this  estimate  are  about  the  same. 
However,  when  we  only  have  discretization  error  we  do  have  F(zh)  —  A(uh',Zh)  —■ ►  0 
as  h  — >  0  and  we  can  make  further  use  of  Z}r  to  determine  how  we  should  adaptively 
refine  a  mesh  so  that  we  can  make  \F(zfl)  —  A(uh',Zh)\  less  than  a  desired  tolerance.  We 
consider  now  how  this  can  be  done  to  adapt  a  space  mesh  when  we  just  have  the  quasi¬ 
static  problem  to  consider. 

In  the  quasi-static  case  the  exact  displacement  U_  satisfies 

a(U_;  ip)  =  0  for  all  appropriate  ip  (5.62) 

(see  (4.19),  i.e.  the  F(.)  term  is  0)  and  the  approximate  displacement^  satisfies 

a{uh  :  ip)  =  0  for  all  appropriate  ip  in  V), .  (5.63) 
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The  dual  solution  in  this  case  is  denoted  by  ip  which  we  assume  is  in  a  space  Vyt  which 
contains  the  finite  element  test  space  14 .  The  estimate  is  now 

J{U)  ~  J(uh)  ~  ~a(uh,±h)  =  -a{Hh,±h  ~  ^±h)  (5-64) 

m 

=  a^t_h  ~  n}th)k  (5-65) 

k= 1 


where  7rW  <G  14  is  the  nodal  interpolant  of  ip  G  14  and  where  M  is  the  number  of  el¬ 
ements.  The  element  quantity  ci(uh,iph  —  ir iph)k  gives  an  indication  of  the  error  in  the 
/.  th  element,  A  =  1,  2,  •  •  •  ,  M.  If  the  aim  of  the  computation  is  to  compute  uh  such  that 
|  J(U_)  —  ■Huh)\  <  5  then  this  is  most  economically  just  achieved  if  we  have 

\a(uh,±h  _7r^fe)fcl  ~  T>>  k  =  1,2,  ■  ■  ■  ,M.  (5.66) 


If  the  element  quantity  for  the  A  th  element  is  larger  than  5/ M  then  we  opt  to  refine  the 
A  th  element.  To  determine  by  how  much  we  should  divide  the  A’th  element  we  consider 
what  is  needed  if  we  are  in  the  asymptotic  convergence  range.  If  we  are  using  degree  p 
polynomials  then  dividing  the  element  into  qy  equal  parts  should  decrease  the  quantity 
corresponding  to  the  region  of  the  Ath  element  by  a  factor  of  about  qf;p .  This  suggests  that 
we  select  qk  so  that 


2 p  _  HHh,±h  -  K±h)k 
%  ~  8/M 


(5.67) 


In  fact  it  is  usually  beneficial  to  take  a  slightly  larger  value  than  this  as  otherwise  a  strategy 
of  aiming  to  almost  exactly  get  to  the  required  accuracy  often  results  in  a  new  mesh  which 
which  has  an  error  which  is  slightly  larger  than  5.  If  the  value  obtained  for  qy  selected  by 
the  above  is  large  then  it  is  also  beneficial  to  replace  it  by  some  fixed  value  with  the  knock 
on  effect  that  more  than  one  mesh  refinement  will  be  needed  before  the  accuracy  of  5  is 
obtained. 


The  mesh  adaption  procedure  used  is  as  follows. 


1.  We  solve  on  a  mesh  of  M  elements  and  compute  the  quantities  a(uh,ip  —  Truy)/,, 
A  =  1,2,  •••  ,  M. 


2.  If 


M 


-  ^±h)k 


k= 1 


<  -5 


then  we  stop  the  computation. 

3.  For  A  =  1, 2,  •  •  •  ,M  we  test  as  follows. 


(5.68) 


If 

then  we  set  qk 


\a{uh,±h-TT±h)k\  < 

1  and  goto  the  next  value  of  A. 


(5.69) 
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Otherwise  we  compute 


Qk 


ceil 


1.05 


a(uh,±h  - 


5/M 


1/(2 p) 


If  qk  >  16  then  we  re-set  qi~  =  16. 

We  divide  the  kth  element  into  qr  parts. 


4.  We  replace  M  by  q\  -f  q->  —  ■■■  +  q^f  and  goto  step  1. 


(5.70) 


The  factor  of  1.05  and  the  bound  of  16  used  above  were  arbitrarily  selected.  These 
values  work  reasonably  well  in  the  test  problems  considered  although  other  values  close 
to  these  are  also  likely  to  work  well. 


6  Results:  discretization  error 

In  this  section  we  test  the  procedure  in  the  case  when  the  error  in  estimating  a  Qol  is  just 
due  to  discretization  error.  The  physical  problem  that  we  are  considering  is  the  nonlinear 
problem  of  the  inflation  of  a  hyperelastic  membrane  introduced  in  section  4.  We  consider 
the  following  two  quantities  of  interest: 

Ji{U)  =  ip(U)  =  [  ip  dr,  =  h0rW  -  ^-{r  +  u\)  (-(r  +  ui)u'3  +  (1  +  u'-^uz)  ,  (6.1) 

Jo  5 

which  is  the  total  potential  energy,  and 

MU)  =  ^  /  rMr)d r  (6.2) 

b  Jo 

which  is  a  measure  of  the  thickness  ratio  A3  in  the  vicinity  of  the  pole  at  r  =  0.  We  take 
b  =  1/8  in  the  computations,  r  =  b  is  a  node  on  each  of  the  meshes  used. 

To  be  specific  for  the  computations,  we  use  a  strain  energy  function  W  of  the  Mooney- 
Rivlin  form  given  by 

fb  =  ^  (^1  +  ^2  +  ^3  —  3)  +  —  (Ax  2  +  A2  2  +  A3  2  —  3)  (6.3) 

where,  as  before,  Ai  and  A2  are  the  principal  stretch  ratios  tangential  to  the  mid-surface  and 
where  A3  =  l/(Ai  A2)  is  the  stretch  ratio  through  the  thickness.  As  the  results  only  depend 
on  the  ratio  P/ho,  where  P  is  the  applied  pressure  and  ho  is  the  undeformed  thickness,  we 
have  taken  ho  =  1  in  our  computations. 

Before  we  consider  estimating  J\  (U_)  and  J2 (U_)  we  show  graphically  various  features 
of  the  solution.  In  figures  6.1  and  6.2  we  show  the  deformed  membranes  at  a  selection  of 
increasing  pressures.  These  plots  depend  on  U_.  In  figures  6.3  and  6.4  we  show  the  stretch 
ratios  Ai,  A2  and  A3  at  the  selection  of  pressures.  These  plots  depend  on  U/ .  The  plots 
suggest  a  'fairly  smooth'  solution  with  not  too  great  a  variation  in  the  quantities  shown 
over  the  domain  0  <  r  <  1.  This  is  confirmed  in  the  results  as  the  meshes  selected  are 
never  too  far  from  being  uniform.  However,  as  we  demonstrate,  the  theory  still  enables 
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us  to  predict  J(U_)  ~  J(lLh )  accurately  and  to  show  how  to  re-mesh  to  achieve  the  required 
accuracy  in  a  small  number  of  steps. 

In  tables  6.1-6.11  we  show  results  obtained  for  estimating  the  total  potential  energy 
in  the  case  of  P  =  1  (in  tables  6. 1-6.5)  and  in  the  case  of  P  =  3  (in  tables  6.6-6.11).  For  each 
pressure  P  we  consider  a  finite  element  space  V),  corresponding  to  linears,  quadratics  and 
cubics  and  in  most  cases  we  consider  the  effect  of  both  uniformly  refining  and  adaptively 
refining.  In  all  the  tables  to  denotes  the  number  of  elements  and  —a(uh.  ph )  is  the  error 
estimate.  As  described  in  section  5.1.3,  the  space  1 4  involves  piecewise  polynomials  of 
one  degree  higher  than  that  which  is  used  for  V),.  In  all  cases  the  so-called  'exact-error' 
is  obtained  by  using  the  finite  element  solution  with  1024  quadratic  element  as  effectively 
the  exact  solution. 

In  all  the  tables,  the  ratio  of  the  error  estimate  and  the  exact  error  are  very  close  to  1 
which  is  consistent  with  the  asymptotic  exactness  described  in  section  3.7.  The  adaptive 
results  described  in  tables  6.2  and  6.4  (for  P  =  1)  and  in  tables  6.7,  6.9  and  6.11  (for  P  = 
3)  indicate  how  well  the  error  estimation  and  adapting  procedure  works  to  achieve  the 
desired  accuracy  in  a  small  number  of  steps.  In  all  cases,  even  for  the  coarsest  mesh  of  8 
elements,  we  appear  to  be  within  the  region  where  the  asymptotic  rate  of  convergence  is 
in  effect.  The  smoothness  of  the  solution  is  evident  by  how  rapidly  the  approximations 
converge  when  the  quadratic  and  cubic  elements  are  used. 

A  similar  batch  of  results  is  obtained  when  J(U_)  =  A  (U_)  and  we  show  these  in 
tables  6.12-6.17. 

In  all  the  above  cases  the  adaptive  procedure  selects  a  mesh  which  is  fairly  uniform. 
As  a  final  test  we  consider  the  case  of  a  non-homogeneous  material  corresponding  to  the 
strain  energy  function 

IF  =  C(r)  j-  (Ai  +  +  A|  -  3)  +  —  (Ax  2  +  A2  2  +  A3  2  -  3)  j  (6.4) 

with 

C(r )  =  1  —  0.9r,  0  <  r  <  1.  (6.5) 

This  describes  a  Mooney  Rivlin  material  which  is  stiffer  at  the  centre  (where  (7(0)  =  1)  than 
at  the  edge  (where  C(  1)  =  0.1).  The  profiles  and  material  paths  in  this  case  for  pressures 
up  to  P  =  1  are  shown  in  figure  6.5.  There  is  now  much  more  tangential  stretching  at 
the  edge  corresponding  to  r  =  1.  With  quadratic  elements  and  a  Qol  which  is  the  total 
potential  energy  we  show  in  figure  6.6  the  variation  of  the  mesh  size  (on  a  log  scale)  with 
r,  0  <  r  <  1,  that  the  procedure  generates  to  achieve  an  accuracy  of  10“ 10  starting  from 
a  uniform  mesh  of  8  elements.  Observe  that  in  the  final  mesh  of  299  quadratic  elements 
the  elements  near  r  =  1  are  about  l/20th  of  the  width  of  those  near  to  the  centre  at  r  =  0. 
Table  6.18  shows  the  error  estimates  generated  by  the  adaptive  procedure  and  this  should 
be  compared  with  table  6.19  which  shows  what  happens  when  we  just  perform  a  sequence 
of  uniform  refinements. 
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Table  6.1:  In  this  table  P  =  1,  we  are  estimating  J\  (U_)  =  —0.050397877719,  we  use  linear 
elements  and  we  perform  uniform  refinement.  The  error  behaves  like  0(h2). 


Elements 

m 

Asymptotic  exactness 

-a(uh,±h)/(Ji(U_)  ~  Ji(uh ))  -  1 

8 

-2.69  x  10"4 

-4.4  x  10”3 

16 

-6.78  x  10"5 

-1.1  x  10"3 

32 

-1.70  x  10-5 

-2.8  x  10~4 

64 

-4.25  x  10"6 

-6.9  x  10“5 

128 

-1.06  x  10-6 

-1.7  x  10“5 

256 

-2.66  x  10-7 

-4.3  x  10"6 

512 

-6.64  x  10-8 

-1.1  x  10“6 

Table  6.2:  In  this  table  P  =  1,  we  are  estimating  J\  (IP),  we  use  linear  elements  and  we 
perform  adaptive  refinement.  The  accuracy  desired  corresponds  to  6  =  10~ ' . 


Elements 

Error  est. 

Asymptotic  exactness 

m 

~a(uhA’h) 

-a(uhAh)/(Ji(U)  ~  JiilLh ))  -  1 

8 

-2.69  x  10-4 

-4.4  x  10"3 

127 

-1.06  x  10-6 

-1.7  x  10"5 

469 

-6.68  x  10"8 

-8.6  x  10"7 

Table  6.3:  In  this  table  P  =  1,  we  are  estimating  J\  (IP),  we  use  quadratic  elements  and  we 
perform  uniform  refinement.  The  error  behaves  like  O ( hr ) . 


Elements 

m 

Asymptotic  exactness 

-a(uh,i’h)/(Ji(U)  ~  JiilLh ))  -  1 

8 

-1.74  x  10"7 

-4.0  x  10"4 

16 

-1.10  x  10"8 

-1.0  x  10"4 

32 

-6.88  x  10-10 

-2.5  x  10-5 

64 

-4.30  x  10~n 

+9.2  x  10"6 

Table  6.4:  In  this  table  P  =  1,  we  are  estimating  .J\  (U_),  we  use  quadratic  elements  and  we 
perform  adaptive  refinement.  The  accuracy  desired  corresponds  to  S  =  10-10. 


Elements 

Asymptotic  exactness 

m 

~a(uhAh)/(Ji(U)  -  Ji(uh))  -  1 

8 

-1.74  x  10~7 

-4.0  x  10"4 

54 

-6.09  x  10"11 

+2.3  x  10"6 

Table  6.5:  In  this  table  P  =  1,  we  are  estimating  Ji(U_),  we  use  cubic  elements  and  we 
perform  uniform  refinement. 


Elements 

m 

Asymptotic  exactness 

~a(uhA’h)/(Ji(U)  -  Ji(uh))  -  1 

8 

-2.69  x  10"11 

-9.0  x  10"4 
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Table  6.6:  In  this  table  P  =  3,  we  are  estimating  J\  (U_)  =  —1.501627519302,  we  use  linear 
elements  and  we  perform  uniform  refinement.  The  error  behaves  like  0(h2). 


Elements 

m 

Error  est. 

~a(uhAh) 

Asymptotic  exactness 

-a{uh,±h)/{Ji{U_)  -  Ji(uh))  -  1 

8 

-7.9  x  10"2 

16 

-1.7  x  10“2 

32 

-6.68  x  10"3 

-4.2  x  10”3 

64 

-1.68  x  10-3 

-1.0  x  10”3 

128 

-4.21  x  10"4 

-2.6  x  10”4 

256 

-1.05  x  10-4 

-6.4  x  10"5 

512 

-2.63  x  10-5 

-1.6  x  10~5 

1024 

-6.58  x  10-6 

-3.8  x  10~6 

2048 

-1.65  x  10-6 

-2.4  x  10”7 

4096 

-4.12  x  10-7 

+2.8  x  10"6 

8192 

-1.03  x  10-7 

+1.2  x  10"5 

16384 

-2.57  x  10-8 

+4.9  x  nr5 

Table  6.7:  In  this  table  P  =  3,  we  are  estimating  .J\  (U_),  we  use  linear  elements  and  we 
perform  adaptive  refinement.  The  accuracy  desired  corresponds  to  5  =  10~ ' . 


Elements 

m 

Asymptotic  exactness 

8 

-7.9  x  10”2 

128 

-4.21  x  10"4 

-2.6  x  10“4 

2030 

-1.65  x  10-6 

-2.4  x  10”7 

9293 

-7.02  x  10"8 

+1.8  x  10~5 
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Table  6.8:  In  this  table  P  =  3,  we  are  estimating  .J\  (U_),  we  use  quadratic  elements  and  we 
perform  uniform  refinement.  The  error  behaves  like  O ( hr ) . 


Elements 

m 

Error  est. 

~a(uh,i>h) 

Asymptotic  exactness 

-a(uh,i’h)/(Ji(U)  -  Jiiuh))  -  1 

8 

-3.29  x  10-4 

-4.05  x  10~3 

16 

-2.09  x  10-5 

-9.86  x  10-4 

32 

-1.31  x  10-6 

-2.44  x  10~4 

64 

-8.22  x  10-8 

-4.58  x  10"5 

128 

-5.14  x  10-9 

+2.30  x  10~4 

256 

-3.21  x  10"1U 

+3.92  x  10~3 

512 

-2.01  x  10-n 

+6.68  x  10-2 

Table  6.9:  In  this  table  P  =  3,  we  are  estimating  J\  (U_),  we  use  quadratic  elements  and  we 
perform  adaptive  refinement.  The  accuracy  desired  corresponds  to  6  =  Id-10. 


Elements 

Error  est. 

Asymptotic  exactness 

m 

-a{uhA’h) 

-a{uh,^h)/{Ji(U)  ~  Ji(uh))  ~  1 

8 

-3.29  x  10-4 

-4.0  x  10-3 

128 

-5.14  x  10-9 

+2.3  x  10-4 

396 

-4.06  x  10_n 

+3.1  x  10-2 

Table  6.10:  In  this  table  P  =  3,  we  are  estimating  we  use  cubic  elements  and  we 

perform  uniform  refinement.  The  error  is  consistent  with  0(h6). 


Elements 

m 

Error  est. 

~a(uh,±h) 

Asymptotic  exactness 

-a(uh,±h)/(Ji(U.)  ~  JiilLh ))  -  1 

8 

-1.27  x  10"6 

-4.6  x  10"3 

16 

-2.03  x  10-8 

-1.1  x  10-3 

32 

-3.20  x  10"1U 

+3.6  x  10-3 

64 

-5.01  x  10"12 

+3.3  x  10_1 

Table  6.11:  In  this  table  P  =  3,  we  are  estimating  we  use  cubic  elements  and  we 

perform  adaptive  refinement.  The  accuracy  desired  corresponds  to  6  =  10-10. 


Elements 

m 

Error  est. 

~a{uhAh) 

Asymptotic  exactness 

~a(uh,^h)/(Ji(U)  ~  Ji(uh))  -  1 

8 

-1.27  x  10"6 

-4.6  x  10"3 

37 

-4.10  x  10"11 

+3.2  x  10"2 
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Table  6.12:  In  this  table  P  =  3,  we  are  estimating  J2(U_),  we  use  linear  elements  and  we 
perform  uniform  refinement.  The  error  behaves  like  0(h2). 


Elements 

m 

Error  est. 

~a(uh,i>h) 

Asymptotic  exactness 

-a{uh,^h)/(J2{U_)  -  J2 (m))  -  1 

8 

-1.46  x  10-3 

-1.0  x  10'1 

16 

-4.12  x  10-4 

-2.1  x  10~2 

32 

-1.06  x  10-4 

-5.1  x  10"3 

64 

-2.67  x  10-5 

-1.3  x  10^3 

128 

-6.70  x  10-6 

-3.2  x  10~4 

256 

-1.68  x  10-6 

-7.9  x  10”5 

512 

-1.9  x  10”5 

1024 

-4.1  x  10"6 

2048 

-2.62  x  10-8 

+2.2  x  10”6 

Table  6.13:  In  this  table  P  =  3,  we  are  estimating  J2(U_)  =  0.046930267582,  we  use  linear 
elements  and  we  perform  adaptive  refinement.  The  accuracy  desired  corresponds  to  S  = 

IQ"7. 


Elements 

Error  est. 

Asymptotic  exactness 

m 

~a{uhAh) 

-a{.uh^h)/{J2(U_)  -  Muh))  ~  1 

8 

-1.46  x  10"3 

-1.0  x  10”1 

128 

-6.70  x  10-6 

-3.2  x  10"4 

1289 

-4.52  x  10-8 

-7.1  x  nr5 
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Table  6.14:  In  this  table  P  =  3,  we  are  estimating  -hiLL),  we  use  quadratic  elements  and  we 
perform  uniform  refinement.  The  error  behaves  like  0(h4). 


Elements 

m 

Asymptotic  exactness 

~a{uh^h)/{MlL)  ~  MUh))  ~  1 

8 

-2.11  x  10-5 

-2.3  x  10-3 

16 

-5.2  x  10-4 

32 

-8.82  x  10-8 

-1.2  x  10-4 

64 

-5.54  x  10-9 

-1.5  x  10-5 

128 

-3.47  x  10"1U 

+2.5  x  10"4 

256 

-2.17  x  10-11 

+4.1  x  10-3 

Table  6.15:  In  this  table  P  =  3,  we  are  estimating  J2(t/),  we  use  quadratic  elements  and  we 
perform  adaptive  refinement.  The  accuracy  desired  corresponds  to  6  =  Id-10. 


Elements 

Asymptotic  exactness 

m 

-a(uh,i’h)/{MU.)  -  Muh))  - 1 

8 

-2.11  x  10-5 

-2.33  x  10~3 

128 

-3.47  x  10_1° 

+2.5  x  10-4 

257 

-1.79  x  10-n 

+4.9  x  10-3 

Table  6.16:  In  this  table  P  =  3,  we  are  estimating  J2(f7),  we  use  cubic  elements  and  we 
perform  uniform  refinement.  The  error  is  consistent  with  0(h6). 


Elements 

Asymptotic  exactness 

m 

-a{uh^h)/{J2{lf)  -  J2(m))  ~  1 

8 

-4.42  x  10"8 

-6.6  x  10"3 

16 

-6.97  x  10_1° 

-1.8  x  10-3 

32 

-1.09  x  10"11 

+7.7  x  10"3 

Table  6.17:  In  this  table  P  =  3,  we  are  estimating  J2(t/),  we  use  cubic  elements  and  we 
perform  adaptive  refinement.  The  accuracy  desired  corresponds  to  5  =  Id-10. 


Elements 

m 

Asymptotic  exactness 

-a(uhA’h)/(J2(U)  ~  Muh))  -  1 

8 

-4.42  x  10"8 

-6.6  x  10"3 

24 

-1.21  x  10'11 

+8.1  x  10"3 
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Figure  6.5:  The  deformed  profiles  and  the  material  paths  for  the  non-homogeneous  sheet 
inflation. 


Table  6.18:  Adaptive  refinements  for  the  non-homogeneous  sheet  inflation  using  quadratic 
elements. 


Elements 

Error  est. 

8 

-1.15  x  10-3 

106 

-2.84  x  10-8 

299 

-3.87  x  10_n 

Table  6.19:  Uniform  refinements  for  the  non-homogeneous  sheet  inflation  using  quadratic 
elements.  _ 


Elements 

Error  est. 

256 

-1.78  x  10“9 

512 

-1.11  x  10"10 

1024 

-6.94  x  10"12 
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7  Results:  modelling  error 


The  graphs  already  shown  in  figure  6.1  give  the  deformed  sheet  profiles  for  different  val¬ 
ues  of  P/ho  assuming  that  the  hyperelastic  membrane  is  deforming  under  quasi-static  con¬ 
ditions,  i.e.  the  deformed  shape  depends  only  on  the  current  load  and  not  on  the  loading 
history  or  the  current  rate  at  which  the  load  is  being  applied.  In  particular  the  quasi-static 
assumption  means  that  inertia  terms  can  be  ignored  which  is  likely  to  be  a  valid  approx¬ 
imation  provided  the  inflation  is  not  done  too  rapidly.  An  obvious  modelling  question  is 
then  'how  slow  does  the  inflation  process  need  to  be  for  a  given  a  quasi-static  solution  to 
be  a  good  approximation  to  a  solution  obtained  for  a  model  in  which  the  inertia  term  is 
included  for  a  given  Qol?'.  A  related  question  is  'how  good  is  the  estimate 


For  the  axisymmetric  inflation  problem  we  are  fortunate  that  it  is  not  too  computationally 
expensive  to  approximately  solve  the  full  dynamic  problem  to  test  the  performance  of  the 
estimate  although  in  other  modelling  situations  this  may  be  an  expensive  computation  if 
a  fine  space  mesh  and  a  large  number  of  time  steps  are  required.  It  must  of  course  not  be 
overlooked  that  to  obtain  the  error  estimate  we  need  to  obtain  ip  and  9h  by  solving  a  dual 
problem  which,  although  linear,  is  a  space  time  problem  which  may  also  need  a  fine  mesh 
and  a  large  number  of  time  steps.  To  summarize  the  steps  involved  we  do  the  following. 


1.  Solve  the  coarse  problem  to  get  uh  and  vh. 

2.  Compute  the  estimate  ). 

3.  Solve  the  dual  problem 


A'( 


4.  Compute  the  error  estimate  F(  (^h^j  )  —  A(  ^  )• 

5.  Test  the  procedure,  if  possible,  by  attempting  to  approximate  U_h  and  V_h  by  approx¬ 
imately  solving  the  fine  problem. 


We  consider  the  same  inflation  problem  as  in  section  6  with  the  same  material.  The 
undeformed  radius  of  the  sheet  is  again  1,  the  initial  thickness  is  taken  as/io  =  10-2  =  0.01 
and  the  density  of  the  incompressible  material  is  taken  as  p  =  1.  Before  the  inflation  is 
started  we  prestretch  the  sheet  so  that  Ci(l)  =  0.1,  the  sheet  is  then  clamped  at  its  edge 
and  then  a  time  dependent  pressure  P(t )  is  applied  where 

P(t)  =  0  <  t  <  T.  (7.2) 

The  final  pressure  is  thus  P(T )  =  0.03  so  that  P(T) /ho  =  3  which  corresponds  to  the 
largest  forcing  action  used  in  figure  6.1  in  the  quasi-static  inflation  case.  For  the  dynamic 
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problem  we  also  need  the  initial  velocity  and  this  is  set  here  as  the  velocity  that  the  quasi¬ 
static  solution  gives  at  t  =  0  so  that  there  is  no  difference  between  the  quasi-static  and 
dynamic  cases  at  this  stage. 

In  the  following  we  consider  the  deformation  that  is  obtained  at  different  final  times 
T  corresponding  to  the  different  values  for  the  rate 


r  = 


0.03 

~T~ 


(7.3) 


for  two  different  Qol.  The  rates  that  are  considered  are  as  follows:  r  =  10~  2 ,  r  =  2  x  10  2 , 
r  =  3  x  10~2,  •  •  ■ ,  r  =  9  x  10-2,  r  =  10  x  10^2,  r  =  20  x  10-2,  r  =  40  x  10~2  and 
r  =  80  x  10-2.  The  final  deformation  at  these  rates  are  shown  in  figure  7.1  and  they  show 
that  there  is  a  large  difference  between  the  profiles  at  the  rates  considered  and  hence  there 
is  a  large  difference  between  the  final  thickness  (i.e.  the  thickness  at  time  t  =  T)  in  these 
cases.  We  consider  the  following  two  quantities  of  interest  concerned  respectively  with 
thickness  and  with  kinetic  energy. 


Dynamic  deformed  profile  at  various  rates  of  pressure  loading 


Figure  7.1:  Deformation  at  final  time  t  =  T  when  the  pressure  is  0.03  corresponding  to  the 
rates  r  =  10“2,  r  =  2  x  10~2,  r  =  3  x  10~2,  ■  ■  ■  ,r  =  9  x  10~2,  r  =  10  x  10-2,  r  =  20  x  10~2, 
r  =  40  x  10”2  and  r  =  80  x  10-2.  The  outermost  curve  is  the  slowest  rate,  the  inner 
parabolic  shaped  curve  is  the  fastest  rate. 
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An  average  thickness  near  the  pole  at  t  =  T 

We  take 


■h(U_)  =  ^  [  r\s(r,T)dr,  6=  f  rdr ,  b  =  0.25  in  the  example,  (7.4) 

o  Jo  Jo 

where  recall  that 

A3  =  — \{  =  {l  +  U[f  +  Ug  and  A2  =  l  +  — • 

A1A2  r 

The  Gateaux  derivative  is  straightforward  to  describe  in  this  case  in  terms  of  quantities 
already  given  in  this  report.  If  we  let  Ai(s),  A 2(s)  and  A 3(s)  denote  the  stretch  ratios  corre¬ 
sponding  to  the  displacement  field  U_  +  sa  then  because  A3  =  1/(AiA2)  we  have,  using  the 
relations  in  equation  (4.41), 


A  '3(U;g) 


-1 


s=0 


-Aq 


Ai—  +  A2 
r 


(AiA2)' 

1  +  U[ 


Ai 


dA2 


ds 


s=0 


+ 


Ai 


/  ^3  / 

al  +  T-  a3 

M 


Thus 


with 


rb 

=  (a  ■  Ja(t,T )  +  a/  •  Ja> (f, T1)]  i'dr 

Jo 


±a(r,T) 


H  0  J’ 


T) 


a|a2  (1  +  Uj\ 
AA,  V  Ui  ) 


t=T 


(7.5) 

(7.6) 


(7.7) 


(7.8) 


A  kinetic  energy  type  term  close  to  the  pole  and  near  the  final  time 


Specifically  we  take 

1  l"T  fb 

■h(V_)  =  -  /  /  r||j/||2  drdf,  with  b  =  0.25. 

2  J o.9T  Jo 

In  this  case  the  Gateaux  derivative  is  simply 


Ji(V_-(3)=  I  I  (3  ■  J_p(r,  t)  rdrdt,  where  J_p(r,  t) 

~  Jo.9T  Jo  ~ 


(7.9) 


(7.10) 


Numerical  results 

In  table  7.1  we  show  the  predictions  for  J3  (U_)  corresponding  to  the  different  rates  at  which 
the  pressure  is  applied.  In  this  case  the  coarse  solution  prediction  (i.e.  the  quasi-static 
solution  prediction)  is  the  same  in  each  case  and  hence  the  error  is  the  difference  between 
each  value  and  0.0478.  The  error  prediction  is  that  obtained  using  the  solution  of  the  dual 
problem.  When  the  rate  is  low  the  prediction  and  the  actual  error  agree  quite  well  (see 
the  rate  1  x  10“  2)  as  the  theory  predicts  whilst  for  the  faster  rates  the  prediction  gives 
the  correct  order  of  magnitude  but  it  is  not  a  sharp  estimate  in  this  case.  This  is  about  as 
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much  as  can  be  expected  as  the  theory  only  says  that  the  estimate  is  good  when  the  coarse 
solution  is  close  to  the  fine  problem.  To  get  the  error  estimate  we  have  to  compute  the 
dual  solution  which  we  find  in  this  example  is  highly  oscillatory  in  both  space  and  time. 
Graphs  of  the  profiles  of  ip i  and  ips  at  the  times  T,  2T/3,  T/3  and  0  and  corresponding 
to  the  fastest  inflation  rate  are  shown  in  figures  7.2,  7.3,  7.4  and  7.5.  The  final  time  value 
of  ip(r,  T)  =  0  in  this  case  but  we  have  0(r.  T)  /  0  because  J3  involves  final  time  values. 
Similar  profiles  are  obtained  with  all  the  other  rates. 


Table  7.1:  The  estimate  of  J3(t/)  by  solving  the  fine  problem,  the  error  by  predicting  using 
the  coarse  problem  and  the  estimate  of  the  error  computed  using  the  dual  solution. 


Rate 

xlO-2 

Estimate 

Error 

Error 

Prediction 

0 

4.782x  10~2 

- 

- 

1 

4.767x  10-2 

-1.477x  10~4 

-1.406xl0”4 

2 

4.782x  10”2 

7.329  xl0“6 

-1.956xl0~4 

3 

5.323  xlO"2 

5.409  xl0“3 

3.647xl0-3 

4 

6.070xl0”2 

1.288xl0-2 

8.I68XIO-3 

5 

7.405  xlO"2 

2.623xl0-2 

1.338xl0-2 

10 

1.438xl0_1 

9.605xl0-2 

3.127x  10-2 

20 

1.775X10”1 

1.297x  10_1 

4.217x  10-2 

40 

2.270x  10-1 

1.792x  10”1 

4.819xl0-2 

80 

4.399  xlO'1 

3.921  xlO"1 

6.511  xl0“2 

In  table  7.2  we  show  the  predictions  for  J^R)  obtained  by  solving  the  fine  prob¬ 
lem,  obtained  by  just  solving  the  coarse  problem  (i.e.  solving  the  quasi-static  problem  for 
the  displacement  and  then  recovering  the  velocity)  together  with  the  actual  error  and  the 
prediction  of  the  error.  In  this  case  the  error  prediction  is  quite  good  throughout  even 
when  the  error  is  large.  The  good  prediction  throughout  seems  to  be  connected  to  a  dual 
solution  which  is  'fairly  simple'  in  that  it  is  not  highly  oscillatory  in  space  or  time.  In  fig¬ 
ures  7.6,  7.7,  7.8  and  7.9  we  show  graphs  of  the  profiles  of  ip  1  and  -03  at  the  times  T,  2 T /3, 
T/3  and  0  respectively.  In  this  case  as  J4 (V_)  does  not  involve  final  time  values  we  have 
ip(r ,  T)  =  9(r ,  T)  =  0,  0  <  r  <  1. 

In  all  the  computations  given  the  coarse  solution  was  obtained  using  M  =  20  quadratic 
elements  and  N  =  200  time  steps.  The  dual  solution  was  obtained  with  40  quadratic  ele¬ 
ments  and  400  time  steps.  This  seems  to  be  enough  (and  possibly  more  than  enough)  to 
ensure  that  the  numbers  given  in  the  tables  here  are  about  accurate  to  all  the  digits  given. 
To  support  this  claim  we  give  in  table  7.3  values  for  the  actual  error  and  the  prediction  of 
this  error  for  the  inflation  rate  of  4  x  10  2  obtained  using  a  different  number  of  space  ele¬ 
ments  and  a  different  number  of  time  steps.  The  closeness  of  the  numbers  in  each  column 
suggests  that  the  space  and  time  meshes  are  fine  enough  and  that  the  error  in  estimating 
the  Qol  in  each  case  is  predominantly  just  the  modelling  error.  The  smoothness  of  the 
deformed  profiles  in  figure  7.1  and  also  the  results  in  table  7.2  indicate  that  not  too  many 
quadratic  elements  are  needed  for  high  accuracy  but  that  a  much  larger  number  of  time 
steps  are  needed  to  get  a  sufficiently  accurate  solution  to  the  fine  problem  and  indeed  to 
the  dual  problem  in  some  cases  (e.g.  the  case  of  predicting  ^3  ( T7 ) ) - 
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Figure  7.2:  Graphs  of  V'l (.,  T)  and  ^3 (.,  T)  at  the  fastest  rate.  Both  are  0. 
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Figure  7.3:  Graphs  of  2T/3)  and  2T/3)  at  the  fastest  rate. 
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Time=1 .650000  zz-thdemo 


Figure  7.4:  Graphs  of  '01  (..  T/3)  and  3)  at  the  fastest  rate. 
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Figure  7.5:  Graphs  of  V’i  (.,  0)  and  -03(.,  0)  at  the  fastest  rate. 
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Table  7.2:  Prediction  of  J4 {V_)  by  solving  the  fine  problem,  by  solving  the  coarse  problem, 
the  error  in  the  prediction  and  the  estimate  of  the  error  by  solving  obtained  using  the 
solution  to  the  dual  problem. 


Rate 

xicr2 

Exact 

Estimate 

Error 

Error 

Prediction 

1 

4.04x  10”2 

4.37x  10”2 

3.32xl0”3 

3-lOxlO”3 

2 

8.07xl0”2 

1.02x10”! 

2.10xl0”2 

1.94xl0”2 

3 

1.21x10”! 

1. 19X10”1 

-2.45  xlO”3 

4.03xl0”3 

4 

1.61xl0_1 

1.02X10”1 

-5.98xl0”2 

-5.33  xlO”2 

5 

2.02x  10”1 

8.12xl0”2 

-1.21x10”! 

-1.23x10”! 

10 

4.03xl0_1 

4.53xl0”2 

-3.58x10”! 

-4.98x10”! 

20 

8.07x10”! 

1.12x  10”2 

-7.96x10”! 

-1.40x10+° 

40 

1.61x10+° 

2.74x  10”2 

-1.59x10+° 

-2.81x10+° 

80 

3.23x10+° 

5.74xl0”2 

-3.17x10+° 

-5.66x10+° 

Table  7.3:  The  'actual  error'  and  the  prediction  of  the  error  using  the  dual  solution  using 
different  numbers  of  quadratic  space  elements  and  time  steps  in  the  case  of  the  inflation 
rate  of  4  x  10”2. 


Thickness  Qol 

K.E.  Qol 

ne 

nt 

Actual  error 

Prediction 

Actual  error 

prediction 

20 

100 

1.284x  10”2 

7.783xl0”3 

-5.960  xlO”2 

-5.440  xlO”2 

20 

200 

1.287x  10”2 

8.165xl0”3 

-5.971  xlO”2 

-5.324  xlO”2 

40 

200 

1.287x  10”2 

8.165xl0”3 

-5.974xl0”2 

-5.326  xlO”2 

40 

400 

1.289xl0”2 

8.169xl0”3 

-5.976  xlO”2 

-5.329  xlO”2 
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Figure  7.6:  Graphs  of  and  ^>3 (.,  T)  at  the  fastest  rate.  Both  are  0. 
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Figure  7.7:  Graphs  of  2T/3)  and  VhG)  2T/3)  at  the  fastest  rate. 
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Figure  7.8:  Graphs  of  '01  (..  T/3)  and  ^{.,T / 3)  at  the  fastest  rate. 
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Figure  7.9:  Graphs  of  V’i  (.,  0)  and  VhO)  0)  at  the  fastest  rate. 
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8  Conclusions 


Computer  simulations  are  an  indispensable  tool  in  modern  engineering  research  and  de¬ 
sign.  They  are  used  to  create  virtual  components  which  can  then  be  subjected  to  virtual 
operating  environments.  The  components  can  be  altered  with  respect  to  geometry,  config¬ 
uration  and  constitution  and  the  simulation  re-virtualised  in  either  the  same  or  a  modified 
virtual  environment.  The  cost  of  simulations  in  silico  is  typically  trivial  compared  with  that 
of  building  hardware  prototypes  and  simulating  their  in-service  actuality  in  a  laboratory. 

The  main  drawback  with  virtual  experiments  is  the  issue  of  reliability:  the  uncer¬ 
tainty  associated  with  how  closely  the  output  matches  that  which  would  be  obtained  in 
the  real  world.  When  simulating  in  the  virtual  world  the  engineer  usually  has  in  mind 
some  Quantity  of  Interest  (Qol),  and  it  is  the  reliable  and  accurate  determination  of  this 
quantity  that  forms  the  goal  of  the  computation.  There  are  two  principal  reasons  why  the 
virtually  realised  Qol  might  not  match  that  which  would  be  observed  in  the  real  world. 

Modelling  error:  this  occurs  when  the  physics  of  a  very  difficult  underlying  problem  are 
approximated  by  those  of  a  simpler  problem  (e.g.  the  approximation  of  a  highly 
heterogeneous  3D  body  by  a  homogenised  2D  cross  section).  It  is  common  to  talk  of 
approximating  a  fine  problem  by  a  coarse  problem,  but  one  should  not  be  led  into 
believing  that  the  coarse  problem  is  necessarily  in  any  sense  easy.  Indeed,  the  entire 
field  of  mathematical  modelling  is,  by  its  nature,  laden  with  modelling  errors  and 
yet  still  contains  many  forbidding  problems. 

Discretization  error:  the  more  manageable  coarser  problem  can  only  be  virtualised  with 
the  aid  of  discretization  techniques  (e.g.  finite  element  methods).  These  generate 
discretization  errors  which  can  be  controlled  and  diminished  by  suitable  enrichment 
techniques  (e.g.  mesh  refinement). 

Methods  for  quantifying  and  controlling  the  second  type  of  error  have  been  studied 
by  numerical  analysts  for  many  years  and  although  for  a  wide  class  of  problems  and  error 
measures,  they  can  be  said  to  be  well  understood,  the  estimation  of  error  in  Qol's  is  a 
young  subject. 

On  the  other  hand,  methods  of  quantifying  and  controlling  the  first  type  of  error  are 
only  partially  understood  for  a  small  class  of  problems.  It  is  a  major  challenge  in  modern 
applied  computational  mathematics  to  devise  theories,  algorithms  and  software  that,  for 
any  given  virtualisation,  can  provide  a  posteriori  a  quantitative  assessment  of  the  size  of 
these  errors  relative  to  the  engineer's  Qol,  and  possibly  adapt  the  virtual  world  in  an 
intelligent  and  automatic  way  so  as  to  reduce  these  errors  to  an  acceptable  level. 

In  this  project  we  have  applied  state-of-the-art  mathematical  theory  to  a  practical  en¬ 
gineering  problem  in  membrane  inflation.  We  have  applied  modern  methods  of  a  posteriori 
error  estimation  in  order  to  address  the  error  in  engineering  quantities  of  interest  that  arise 
from  modelling  simplifications  and  also  from  finite  element  approximation.  Our  results 
show  that  this  technique  works  in  that  it  can  be  used  to  adaptively  optimise  the  finite  ele¬ 
ment  mesh,  as  well  as  give  an  indication  of  the  physical  error  contained  in  the  modelling 
assumptions. 
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The  eventual  benefit  of  research  of  this  nature  and  the  downstream  uptake  of  the  tech¬ 
niques  and  methodology  into  commercial  engineering  software  cannot  be  overstated.  By 
isolating  and  estimating  the  two  principal  sources  of  error  in  physical  simulations  (mod¬ 
elling  and  discretization)  the  software  will  allow  the  practitioner  to  assign  a  degree  of 
confidence  in  the  results  that  has  never  before  been  possible. 

However,  what  is  absolutely  clear  from  our  efforts  thus  far  is  that  to  make  further 
progress  in  research  of  this  nature  some  form  of  'computational  workhorse'  is  needed.  So 
far,  as  is  common  in  primary  research,  we  have  generated  our  own  software  using,  where 
possible,  codes  from  earlier  efforts.  This  approach  has  certain  advantages  in  terms  of  code 
control  and  ownership,  but  has  severe  disadvantages  in  terms  of  the  time  it  takes  to  bring 
a  theoretical  result  to  practical  fruition. 
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9  Appendix:  Some  of  the  notation 


Xi,  x2,  x3 

Xi,  X2,  X3 

r,  ip,  X3 

§li>  Xi,  e3, 

&r'  &ip' 

t 

^3  D 

n 

T 

Q 

i 

n 

u 

u 

V 

v 

V 


Cartesian  coordinates.  Typically  upper  case  letter  are  used  when  we  con¬ 
sider  the  undeformed  position  of  a  point. 

Cartesian  coordinates.  Typically  lower  case  letter  are  used  when  we  con¬ 
sider  the  deformed  position  of  a  point. 

Cylindrical  polar  coordinates,  r  is  the  radius,  ip  is  the  angle  and  X3  is  the 
vertical  coordinate. 

These  denote  the  cartesian  base  vectors  when  cartesian  coordinates  are 
being  considered. 

These  denote  the  base  vectors  when  cylindrical  polars  are  being  consid¬ 
ered. 

t  denotes  time  when  space-time  problems  are  considered. 

The  region  of  a  general  three-dimensional  body  when  it  is  undeformed. 
The  region  of  the  undeformed  mid-surface  of  the  membrane. 

The  final  time. 

This  is  used  for  the  space  time  region  Q  =  tt  x  ( 0,  T). 

When  a  function  only  depends  on  one  variable  this  is  used  to  denote  the 
derivative  with  respect  to  that  variable.  This  is  extensively  used  in  axi- 
symmetric  problems  when  we  have  quantities  which  only  depend  on  r. 
The  second  derivative  of  a  function  of  one  variable.  Similarly  denotes 
the  third  derivative  etc.. 

The  first  time  derivative. 

The  second  time  derivatives. 

The  exact  displacement.  When  cartesian  coordinates  are  being  used  U  = 
(LJ\ ,  U2,  U3)t  with  Ui  denoting  the  component  in  the  e,  direction. 

When  cylindrical  polar  coordinates  are  being  used  and  axi-symmetric  de¬ 
formations  are  being  considered  we  let  U_  =  (U\ .  U3)r  with  U\  denoting 
the  component  in  the  er  direction  and  with  U3  denoting  the  component  in 
the  direction. 

An  approximate  displacement.  When  cartesian  coordinates  are  being  used 
u  =  (ui,u2,u3)t  with  u,  denoting  the  component  in  the  e,  direction. 
When  cylindrical  polar  coordinates  are  being  used  and  axi-symmetric  de¬ 
formations  are  being  considered  we  let  u  =  (u±,u3)T  with  u\  denoting  the 
component  in  the  er  direction  and  with  u3  denoting  the  component  in  the 
e3  direction. 

The  exact  velocity.  In  the  application  we  enforce,  weakly,  the  condition 
V  =  U_.  When  cylindrical  polar  coordinates  are  being  used  and  axi- 
symmetric  deformations  are  being  considered  we  let  V_  =  (VI ,  V3)r  with 
W\  denoting  the  component  in  the  er  direction  and  with  W3  denoting  the 
component  in  the  e3  direction. 

An  approximate  velocity.  When  cylindrical  polar  coordinates  are  being 
used  and  axi-symmetric  deformations  are  being  considered  we  let  v  = 
(,vi,v3)t  with  v\  denoting  the  component  in  the  er  direction  and  with  v3 
denoting  the  component  in  the  e3  direction. 

In  the  abstract  formulation  V  denotes  a  function  space. 
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J(.) 

p 

ho 

P 

V3  d 


V 


Div 


l3D,I 


F3  D 

F 

Ai,  A2,  A3 


W 

a 

II3D 


The  quantity  of  interest  (Qol).  This  will  typically  depend  on  the  displace¬ 
ment  and/ or  the  velocity  and/or  the  first  spatial  derivatives  of  these. 

The  density  of  the  material. 

The  undeformed  thickness  of  a  sheet. 

The  applied  pressure.  Typically  we  have  P  =  P(t). 

The  gradient  operator  in  the  fully  three-dimensional  case. 


v3Du 


dU  du  dU  \ 
dJ^’dx^’dx^J ' 


The  gradient  operator  in  the  membrane  case.  When  the  membrane  model 
is  being  considered  VC/  means  the  3x2  tensor 


W 


dU  dU  \ 
dXl'dX^)' 


The  divergence  operator  with  respect  to  X\,  X2  and  X:> .  In  the  context 
when  equations  of  motion  and/or  quasi-static  equilibrium  are  being  de¬ 
scribed  and  we  have  the  fully  three  dimensional  case,  the  divergence  of 
U3D  is  the  vector 


Div  n3D 


(diin 

+ 

dn21 

dX2 

+ 

dn.{i  \ 
cDC 

du12 

+ 

an22 

+ 

dn32 

sx^ 

an  13 

+ 

sn23 

+ 

dn33 

uxx2 

Identity  tensors 


1 3D 


The  3x3  deformation  gradient  when  a  fully  three  dimensional  defor¬ 
mation  is  being  considered.  It  is  the  Jacobian  matrix  of  the  mapping 
X  — >  x  =  X_  +  U_  and  is  given  by  F3D  =  I3D  +  VC/. 

The  3  x  2  membrane  deformation  gradient.  In  the  membrane  model  ap¬ 
plication  this  is  the  first  two  columns  of  F^d. 

The  principal  stretch  ratios,  i.e.  A?,  i  =  1,2,3,  are  the  eigenvalues  of 
F^DF3£>  and  of  FsdFJd  in  the  fully  three-dimensional  case.  In  the  mem¬ 
brane  case  A j,  i  =  1,2,  are  the  eigenvalues  of  FTF  and  correspond  to 
stretching  in  directions  tangential  to  the  mid-surface.  In  the  incompress¬ 
ible  membrane  case  we  have  A3  =  1/ (Ai A2)- 

The  strain  energy  function  of  a  hyperelastic  material.  In  the  case  of  incom¬ 
pressible  membrane  models  we  have  W  =  W(Ai,  A2). 

The  3  x  3  Cauchy  stress  tensor.  This  is  symmetric. 

The  3  x  3  nominal  stress  tensor.  This  is  not  in  general  symmetric  and  it  is 
related  to  a  by  the  relation  II32}  =  (det  F3d)F3£,0-. 
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n7 

«(■;■) 

«(■;•) 

F(-) 


a' A>  (■.,.) 


F'U-) 


h 1(n) 


The  2  x  3  nominal  stress  tensor  in  the  membrane  model.  As  a  consequence 
of  the  membrane  simplification  this  is  the  first  two  rows  of  II3  o  ■ 

The  3  x  2  first  Piola  stress  tensor  in  the  membrane  model. 

The  double  dot  product  operation,  i.e.  for  tensors  A  and  B  of  compatible 
sizes  we  have  A  :  B  =  •  AijBij  =  tr(A7  B). 

The  form  used  in  the  description  of  the  quasi-static  problem.  The  semi¬ 
colon  indicates  that  the  functional  is  non-linear  in  the  term  before  the  semi¬ 
colon  and  is  linear  in  the  term  afterwards. 

The  form  used  in  the  abstract  description  and  also  the  form  used  for  the 
weak  form  of  the  dynamic  problem. 

In  the  formulations  this  is  used  for  the  integrand  or  part  of  the  integrand 
in  the  expressions  for  a(.; .)  and/or  for  A(.; .). 

A  linear  functional  appearing  on  the  right  hand  side  in  expressions  such 


as  a(U_ ;  ip) 


F{ip )  or  A{ 


ip 


7T  )  =  n  7T 


ip 


where  ip  and  9  are  test 


vectors. 

The  first  Gateaux  derivatives  of  a  and  A  with,  for  example. 


a'(U_;  a,  ip)  = 


ds 


a(U_  +  sa-  ip) 


s= 0 


s=0 


The  term  just  to  the  right  of  the  semi-colon  indicates  the  "direction"  in 
which  we  differentiate.  Higher  derivatives  are  similarly  defined. 

The  first  Gateaux  derivative  of  F.  F'(U_ ;  ip)  means 


F\U]ip)  =  -t-F(U_  +  sip) 
—  as  — 


s= 0 


Higher  derivatives  are  similarly  defined. 

The  unknown  displacement  and  velocity  in  the  dynamic  problem. 


The  test  vector  in  the  dynamic  problem  and  the  unknown  vector  in  the 
dual  problem. 

The  test  vector  in  the  dual  problem. 

The  space  of  functions  which  are  square  integrable  over  the  domain  ft . 
The  space  of  functions  whose  first  weak  derivatives  are  square  integrable 
over  the  domain  fl 

Ui(x,t)  is  in  this  space  if  Ui(.,t)  £  fT1(ff)  for  each  time  t  €  (0,  T)  and 
Ui(x,.)  G  Hl(0,T)  for  each  point  x  £  ft. 

Vi(x,t )  is  in  this  space  if  Vi(.,t)  £  ^(fl)  for  each  time  t  £  (0,  T)  and 
Vi(x, .)  £  Hl{ 0,  T)  for  each  point  x  £  ft. 

ipi(x,  t)  is  in  this  space  if  ipi(.,t)  £  H1(ft)  for  each  time  t  £  (0,  T)  and 
ipi(x, .)  £  1/2(0,  T)  for  each  point  x  £  fl 


-ff1((0,  T),H1(fl)) 
H'd  0,T),L2(fi)) 
£2((o,  T),  H1  (fi)) 
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L2((0,T),L2(Q)) 


6i(x,t )  is  in  this  space  if  9i(.,t)  G  for  each  time  t  G  (0,  T)  and 

8i(x,.)  G  £2(0,  T)  for  each  point  x  G  fi. 
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